% EXAMPLE2_42 Computations for Example 2.42.
%
clear r mc mp iq idelta irout AB
global r mc mp iq idelta irout AB
f1=' %4.0f %18.10e %18.10e %10.2e %9.2e\n';
f2=' %53.2e %9.2e\n';
mc=1; mp=0; iq=1; idelta=2; irout=1; AB=[0 Inf];
N=20; Mmax=300; eps0=1e2*eps;
load -ascii abeinstein1
load -ascii abeinstein2
load -ascii abfermi1
load -ascii abfermi2
abexe1=abeinstein1(1:N,:);
abexe2=abeinstein2(1:N,:);
abexf1=abfermi1(1:N,:);
abexf2=abfermi2(1:N,:);
fprintf('\n')
disp('    n       alpha_n            beta_n           erra      errb')
r=1;
[ab,Mcap,kount]=mcdis(N,eps0,@quadeinstein,Mmax);
errabm=max(abs(ab-abexe1)./moll(abexe1));
n=[5 10 20];
errab=abs((ab(n,:)-abexe1(n,:))./moll(abexe1(n,:)));
abe=[(n-1)' abexe1(n,:) errab];
fprintf(f1,abe')
fprintf(f2,errabm')
r=2;
[ab,Mcap,kount]=mcdis(N,eps0,@quadeinstein,Mmax);
errabm=max(abs(ab-abexe2)./moll(abexe2));
errab=abs((ab(n,:)-abexe2(n,:))./moll(abexe2(n,:)));
abe=[(n-1)' abexe2(n,:) errab];
fprintf(f1,abe')
fprintf(f2,errabm')
fprintf('\n')
r=1;
[ab,Mcap,kount]=mcdis(N,eps0,@quadfermi,Mmax);
errabm=max(abs(ab-abexf1)./moll(abexf1));
errab=abs((ab(n,:)-abexf1(n,:))./moll(abexf1(n,:)));
abe=[(n-1)' abexf1(n,:) errab];
fprintf(f1,abe')
fprintf(f2,errabm')
r=2;
[ab,Mcap,kount]=mcdis(N,eps0,@quadfermi,Mmax);
errabm=max(abs(ab-abexf2)./moll(abexf2));
errab=abs((ab(n,:)-abexf2(n,:))./moll(abexf2(n,:)));
abe=[(n-1)' abexf2(n,:) errab];
fprintf(f1,abe')
fprintf(f2,errabm')
