% INTEXP0INF Integration relative to the weight function w(t)=exp(-1/t-t)
%
f1='%5.0f %21.13e\n';
fprintf('\n')
disp('    n      n-point Gauss')
load -ascii abexp0inf;
ab=abexp0inf; 
for n=4:4:24
  xw=gauss(n,ab);
  int=sum(xw(:,2).*mfun('BesselJ',0,xw(:,1)));
  fprintf(f1,n,int)
end
fprintf('\n')
abl=r_laguerre(1000);
for n=100:100:1000
  xwl=gauss(n,abl);
  intl=sum(xwl(:,2).*mfun('BesselJ',0,xwl(:,1))...
    .*exp(-1./xwl(:,1)));
  fprintf(f1,n,intl)
end 
fprintf('\n')
abl=r_laguerre(40);
load -ascii abexp0;
abe=abexp0;
for n=4:4:40
  xwe=gauss(n,abe); xwl=gauss(n,abl);
  inte=sum(xwe(:,2).*mfun('BesselJ',0,xwe(:,1))...
    .*exp(-xwe(:,1)))
  intl=exp(-1)*sum(xwl(:,2).*mfun('BesselJ',0,...
  1+xwl(:,1)).*exp(-1./(1+xwl(:,1))))
  int=inte+intl;
  fprintf(f1,n,int)
end