%SR_TWOINT

function ab=sr_twoint(dig,N,xi,p,q)
digits(dig); sxi=vpa(xi); sp=vpa(p); sq=vpa(q);
ab=vpa(zeros(N,2)); sxi2=sxi^2; s1=vpa(1); s2=vpa(2);
ab(1,2)=(s1-sxi2)^(sp+sq+s1)*gamma(sp+s1)*gamma(sq+s1)/ ...
  gamma(sp+sq+s2);
if N==1, return; end
sabj=sr_jacobi(dig,N,sq,sp);
ab(2,2)=vpa(1/2)*(s1-sxi2)*sabj(1,1)+vpa(1/2)*(s1+sxi2);
if N==2, return; end
for k=1:floor((N-1)/2)
  ab(2*k+1,2)=(vpa(1/2)*(s1-sxi2))^2*sabj(k+1,2)/ab(2*k,2);
  ab(2*k+2,2)=vpa(1/2)*(s1-sxi2)*sabj(k+1,1)+vpa(1/2)*(s1+sxi2)-ab(2*k+1,2);
end
