% EXAMPLE3_37 Computations for Example 3.37.
%
clear eps0 Nmax M theta
global eps0 Nmax M
global theta
%
% Be sure to declare the variable theta global also in the routines
% R_MOD and QUADRAT and make the changes indicated in the routines 
% QUADRAT and FEX31.
%
f1='            eta=%4.3f theta=%8.2e k=%2.1f eps0=%8.2e\n\n';
f2='%4.0f %4.0f %22.14e %12.4e %5.0f %5.0f %5.0f\n';
fprintf('\n')
disp('                           M=2*N-1')
%disp('                      M=2*floor((N+1)/2)-1')
%disp('                             M=1')
%disp('                             M=0')
Nm=10; eps0=100*eps; Nmax=200; numax=50;
eta=-0.001; theta=1e-4; 
k=.5;
%k=1.5;
%k=2.5;
exact=2.21715010091142;
%exact=1.78001232862939;
%exact=3.74038445836967;
%
% When theta is "large", increase eps0 to avoid an excessively
% large value of Nmax that would be needed for convergence. Also
% increase Nmax appropriately. For example,
%            eps0=1e6*eps; Nmax=800;
%
fprintf(f1,eta,theta,k,eps0)
disp('   n    m        integral            error      Ncap  kount  nu')
for N=1:Nm
  M=2*N-1;
%  M=2*floor((N+1)/2)-1;
%  M=1;
%  M=0;
  [xw,Ncap,kount,nu]=bose_einstein_diffpole(N,M,eta,theta,k,eps0,Nmax,numax);
  int=sum(xw(:,2).*fex37(xw(:,1),eta,theta));
  err=abs((int-exact)/exact);
  if ~all(N-[1 4 7 10])
    fprintf(f2,N,M,int,err,Ncap,kount,nu)
  end
end
