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