%EXAMPLE2
%
f0='%10.0f %8.4f %4.0f %6.0f %18.12f\n';
dig=48; digits(dig);
N=28;
ab0=sr_jacobi01(dig,N,-1/2,0);
disp('       beta     c      m      n           I') 
for beta=[100 200 400]
  for c=[1.01 1.001 1.0001]
    for m=[0 1 3 5 7 9]
      if m==0
        ab=subs(ab0);
      else
        z(1)=vpa(c);
        z(2)=vpa(pi*i/beta); z(3)=conj(z(2));
        z(4)=vpa(3*pi*i/beta); z(5)=conj(z(4));
        z(6)=vpa(5*pi*i/beta); z(7)=conj(z(6));
        z(8)=vpa(7*pi*i/beta); z(9)=conj(z(8));
        ab=mod_md(N,m,z,ab0,@rho0);
        ab=subs(real(ab));
      end
      for n=4:4:N
        xw=gauss(n,ab); x=xw(:,1); w=xw(:,2);
        I=sum(w.*f2(x,m,beta,c));
        fprintf(f0,beta,c,m,n,I)
      end
      fprintf('\n')
    end
  end
 end
