function lindh_integ_max2
%This is the calculation of the Lindh Gaussian - Eq. (3.11) with Maxwell
%quadrature points (p = 2) - Figs. 3.4(A) and 3.4(B) 
format long e
hold off
% Scale factor range
alfa=4; smin=1/sqrt(alfa)-0.1; smax=1/sqrt(alfa)+0.6; ds=0.0002;
scx=1/sqrt(alfa);  ell=0;
% There are two choices for the order of the Maxwell quadrature
% namely, 2 and 8, and 4 and 10. 
%order=[2,8];
order=[4 10]; xm=[]; kmax=0; yerr=[];
for k=1:1:2
       m=order(k);
%       fprintf(myfile,'%i\n',m);
ptwt=maxwell2(m);
x=ptwt(:,1); w=ptwt(:,2); n1=length(x);
       xm=[xm m]; sc=[]; yerr=[];
% Calculate the integral versus the scale factor       
for scale=smin:ds:smax
    sc=[sc scale];
% This is the small loop through N - 3 curves
% This is the larger loop through s
% See Eq. (3.3)
bigw=w./(x.^2.*exp(-x.*x)); ws=scale*bigw; xs=scale*x;
fx=xs.^(ell+2).*exp(-alfa*xs.^2);
%fx=10*xs.^2.*exp(-10*xs.^2);
s=2*alfa^((ell+3)/2)/gamma((ell+3)/2)*sum(ws.*fx);
err=log10(abs(s-1.d0)/1.d0);
%fprintf(myfile,'%16.8f %13.5e\n',scale,s-1);
yerr=[yerr err];
end
kmin=kmax+1; kmax=kmin+1;
plot(sc,yerr,'-k','linewidth',1.6)
hold on
xlabel('${\rm s}$','Interpreter','latex','fontsize',32)
ylabel('$\log_{10}[{\rm Relative}\;\; {\rm Error}]$','Interpreter','latex','fontsize',32)
axis([smin smax -16 0])
set(gca,'FontSize',32)
set(gca,'Ytick',[-16:4:0],'linewidth',1.6)
set(gca,'Xtick',[smin:.2:smax],'linewidth',1.6)
end
% Format for the first graph
%Graph (a)
%str1={'N = 2'}; %str2={'N = 8'}; %str3={'$\alpha = 4$'}
%text(0.55,-5,str1,'Interpreter','latex','fontsize',32)
%text(0.8,-8,str2,'Interpreter','latex','fontsize',32)
% ======================================================
% Format for the second graph
%Graph (b)
str1={'N = 4'}; str2={'N = 10'}; str3={'$\alpha = 4$'};
text(0.66,-6,str1,'Interpreter','latex','fontsize',32)
text(0.7,-10,str2,'Interpreter','latex','fontsize',32)
% =======================================================
text(0.9,-14,str3,'Interpreter','latex','fontsize',32)
set(gcf, 'Units','centimeters','Papersize',[48,48])
set(gcf, 'Units','centimeters','Position',[3 3 24 20])
% Maxwell quadrature points and weights for p = 2
function ptwt = maxwell2(m)
load abmaxp2.dat; n=90;
a=abmaxp2(1:m,1); b=abmaxp2(1:m,2); rtb=sqrt(b); rtb(m)=[];
t=diag(rtb,-1)+diag(a)+diag(rtb,1);
[f,lambda]=eig(t);
pt=diag(lambda);
wt=sqrt(pi)*f(1,:).^2/4.d0;
ptwt=[pt,wt'];