%R_TWOINT

function ab=r_twoint(N,xi,p,q)
ab=zeros(N,2); xi2=xi^2;
ab(1,2)=(1-xi2)^(p+q+1)*factorial(p)*factorial(q)/ ...
  factorial(p+q+1);
if N==1, return; end
abj=r_jacobi(N,q,p);
ab(2,2)=.5*(1-xi2)*abj(1,1)+.5*(1+xi2);
if N==2, return; end
for k=1:floor((N-1)/2)
  ab(2*k+1,2)=(.5*(1-xi2))^2*abj(k+1,2)/ab(2*k,2);
  ab(2*k+2,2)=.5*(1-xi2)*abj(k+1,1)+.5*(1+xi2)-ab(2*k+1,2);
end
