%NEUTRAL_EX4
%
f1='%8.0f %4.0f %20.15f\n';
disp('       n    m            I') 
global imod
for dig=32
%for dig=16
  digits(dig);
  z=-vpa(1/1000); alpha=-1/2;
  N=50; m=0;
  ab0=sr_jacobi01(dig,N,alpha,0);
  for n=10:10:N
    xw=gauss(n,subs(ab0));
    x=xw(:,1); w=xw(:,2);
    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=10
    ab0=sr_jacobi01(dig,N,alpha,0);
    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.*cos(pi*x/2).*(x-subs(z)).^(m+alpha));
      fprintf(f1,n,m,I)
    end
    fprintf('\n') 
  end
end
