% EXAMPLE3_36 Computations for Example 3.36.
%
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 FEX37.
%
f1='            eta=%4.1f theta=%8.2e k=%2.1f eps0=%8.2e\n\n';
f2='%4.0f %4.0f %22.14e %12.4e %7.0f %7.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;
eta=-1; theta=1e-4;
%eta=-1; theta=1;
k=.5;
%k=1.5;
%k=2.5;
exact=.37970886599807399070;
%exact=.52608888707964629059;
%exact=1.2665691265431175469;
%
% The next values of exact relate to theta=1, theta=10,
% and theta=100.
%
%exact=.487530636003353;
%exact=1.00386902661383;
%exact=2.9219477617;
%
% 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')
for N=1:Nm
  M=2*N-1;
%  M=2*floor((N+1)/2)-1;
%  M=1;
%  M=0;
  [xw,Ncap,kount]=bose_einstein(N,M,eta,theta,k,eps0,Nmax);
  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)
  end
end
