function [I,N,can]=sintegral(dig,n,bet,gam,eps0)
digits(dig);
Nmax=100;
c=vpa(sqrt(1+(gam/2)^2));
r0=vpa((2/gam)*sqrt(1+c));
rho0=vpa(2*r0*atan(r0));
ab0=sr_jacobi01(dig,Nmax,-1/2,0);
if gam<=1/2, ab=schri4(dig,Nmax,ab0,c,rho0); end
N=0;
I0=vpa(0); I=vpa(1);
err=abs(I-I0);
while subs(err)>eps0
  N=N+5; 
  if N>Nmax
    error('N exceeds Nmax')
  end
  I0=I;
  if gam>1/2
    xw=sgauss(dig,N,ab0(1:N,:)); x=xw(:,1); w=xw(:,2);
    if n==2*floor(n/2)
      fun=-(1/4)*(w.*(schebpol(dig,x,n)./((x+c).*sqrt(1+x)))./ ...
          (x-c));
    else
      fun=(1/4)*(w.*(sinh(bet*x)./(1+cosh(bet*x))).* ...
          (schebpol(dig,x,n)./((x+c).*sqrt(1+x)))./(x-c));
    end
  else
    xw=sgauss(dig,N,ab(1:N,:)); x=xw(:,1); w=xw(:,2);
    if n==2*floor(n/2)
      fun=-(1/4)*(w.*(schebpol(dig,x,n)./((x+c).*sqrt(1+x))));
    else
      fun=(1/4)*(w.*(sinh(bet*x)./(1+cosh(bet*x))).* ...
          (schebpol(dig,x,n)./((x+c).*sqrt(1+x))));
    end
  end
  I=sum(fun); M=max(abs(subs(fun))); 
  can=log10(M/abs(subs(I)));
  err=abs(I-I0);
  if abs(subs(I))>1, err=err/abs(I); end
end
