%NEUTRAL_EX1
%
f0='%8.0f %4.0f %14.3e\n';
f1='%8.0f %4.0f %20.14f\n';
disp('       n    m        error')
%disp('       n    m            I') 
global imod
for dig=28
  digits(dig);
  z=-vpa(1/1000); alpha=-1/2;
  N=50; m=0;
  ab0=sr_jacobi01(dig,N);
  Iexact=1.937754196921554;
  for n=10:10:N
    xw=gauss(n,subs(ab0));
    x=xw(:,1); w=xw(:,2);
    I=sum(w.*(x-subs(z)).^alpha);
    error=abs(I-Iexact);
    fprintf(f0,n,m,error)
%    I=sum(w.*cos(pi*x/2).*(x-subs(z)).^alpha);
%    fprintf(f1,n,m,I)
  end
  fprintf('\n')
  imod=1;
  for m=24
%  for m=18:-6:6
%  for m=24:6:42
    if m>24, dig=52; end
    ab0=sr_jacobi01(dig,N);
    z=-vpa(1/1000);
    ab=mod_d(N,m,z,ab0,@rho01_mu);
    ab=subs(ab);
    for n=10:10:N
      xw=gauss(n,ab);
      x=xw(:,1); w=xw(:,2);
      I=sum(w.*(x-subs(z)).^(m+alpha));
      error=abs(I-Iexact);
%      I=sum(w.*cos(pi*x/2).*(x-subs(z)).^(m+alpha));
      fprintf(f0,n,m,error)
%      fprintf(f1,n,m,I)
    end
    fprintf('\n')
  end
end
