% 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')