%EXAMPLE3
%
global r
f0='%10.0f %6.0f %18.12f %14.4e\n';
f1='%16.0f %18.12f %14.4e\n';
dig=48; digits(dig);
N=12; 
ab0=sr_laguerre(dig,N);
disp('         m      n           I             error') 
for r=1:2
  for m=[0 2 4]
    if m==0
      ab=subs(ab0);
    else
      z(1)=vpa(2*r*pi*i); z(2)=conj(z(1));
      z(3)=vpa(4*r*pi*i); z(4)=conj(z(3));
      ab=mod_md(N,m,z,ab0,@rho0);
      ab=subs(real(ab));
    end
    for n=2:2:N
      xw=gauss(n,ab); x=xw(:,1); w=xw(:,2);
      I=sum((1/r)*w.*exp(-x/r).*f3(x,m));
      if r==1
        err=abs(I-.6449340668482264);
      elseif r==2
        err=abs(I-.4816405210580757);
      end
      fprintf(f0,m,n,I,err)
    end
    fprintf('\n')
  end
  load -ascii abeps1;
  abe1=abeps1;
  load -ascii abeps2;
  abe2=abeps2;
  for n=2:2:N
    if r==1
      xwe=gauss(n,abe1); xe=xwe(:,1); we=xwe(:,2);
    elseif r==2
      xwe=gauss(n,abe2); xe=xwe(:,1); we=xwe(:,2);
    end
    Ie=sum(we.*exp(-xe)); 
    if r==1
      erre=abs(Ie-.6449340668482264);
    elseif r==2
      erre=abs(Ie-.4816405210580757);
    end
    fprintf(f1,n,Ie,erre)
  end
  fprintf('\n')
end
