% EXAMPLE2_43 Computations for Example 2.43.
%
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=2; mp=0; iq=1; idelta=2; irout=1; AB=[[0 10];[10 Inf]];
N=20; Mmax=300; eps0=1e2*eps;
load -ascii absqp1einstein1;
load -ascii absqp1einstein2;
load -ascii absqp1fermi1;
load -ascii absqp1fermi2;
abexe1=absqp1einstein1(1:N,:);
abexe2=absqp1einstein2(1:N,:);
abexf1=absqp1fermi1(1:N,:);
abexf2=absqp1fermi2(1:N,:);
fprintf('\n')
disp('    n       alpha_n            beta_n           erra      errb')
r=1;
[ab,Mcap,kount]=mcdis(N,eps0,@quadsqeinstein,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,@quadsqeinstein,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,@quadsqfermi,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,@quadsqfermi,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')