% SI101 Symbolic version of I101.m
%
function y=sI101(dig,N,n0,a,b)
eps0=.5*10^(-dig); eta=1/100;
k0=k_0(eps0,a,b); 
[k0]
dig=dig+4;
tic
sab=sr_jacobi01(dig,N);
err=1; ss1=vpa(0,dig);
n=n0;
while err>eps0
  ss0=ss1; n=n+5;
  if n>N
    error('n exceeds N in sI101.m')
    break
  end
  sxw=sgauss(dig,n,sab(1:n,:));
  ss1=vpa(0,dig);
  for k=1:n
    st=vpa(eta+(1-eta)*sxw(k,1),dig);
    sy=vpa(-st*exp(-st),dig);
    sw=swofy1(dig,sy);
    sw=(-sw)^vpa(a,dig)*st^(-vpa(b,dig));
    ss1=vpa(ss1+sxw(k,2)*sw,dig); 
  end
  ss1=vpa((1-eta)*ss1,dig);
  err=subs(abs((ss1-ss0)/ss1));
  [err]
end
[n ss1]'
t1=toc;
tic
int=vpa(0,dig);
for k=2:k0+1
  if k<=k0
    c=vpa(10^(-(k+1)),dig); d=vpa(10^(-k),dig);
  else
    c=vpa(0,dig); d=vpa(10^(-k),dig);
  end
  err=1; int1=vpa(0,dig); 
  n=0; nmax=0;
  while err>eps0
    n=n+5;
    int0=int1;
    sxw=sgauss(dig,n,sab(1:n,:));
    ss=vpa(0,dig);
    for nu=1:n
      sx=vpa(c+sxw(nu,1)*(d-c),dig);
      sw=swofy1(dig,-sx*exp(-sx));
      sw=(-sw)^vpa(a,dig)*sx^(-vpa(b,dig));
      ss=vpa(ss+sxw(nu,2)*sw,dig);
    end
    int1=vpa((d-c)*ss,dig);
    err=abs(subs(int1-int0));
  end
  if n>nmax
    nmax=n;
  end
  int=vpa(int+int1,dig);
  [k nmax int1 int]'
end
y=vpa(ss1+int,dig);
t2=toc;
[t1 t2]'