% EXAMPLE3_35 Computations for Example 3.35.
%
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 FEX35.
%
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')
%disp('                      M=2*floor((N+1)/2)')
%disp('                             M=2')
%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=.29051241701949266261;
%exact=.46087845417799195532;
%exact=1.1860735010757557836;
%exact=1.3964418203491153395;
%exact=2.6618732791071501381;
%exact=7.6272560956534476326;
%
% The next values of exact relate to theta=1, theta=10,
% and theta=100, respectively.
%
%exact=.383869768812140;
%exact=.820188540209559;
%exact=2.41632871284393;
%
% 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;
%  M=2*floor((N+1)/2);
%  M=2;
%  M=0;
  [xw,Ncap,kount]=fermi_dirac(N,M,eta,theta,k,eps0,Nmax);
  int=sum(xw(:,2).*fex35(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