function ab=sr_jacobi(dig,N,a,b)
% SR_JACOBI Symbolic counterpart of r_jacobi.m
%
digits(dig);
if nargin<3, a=0; end;  if nargin<4, b=a; end
if((N<=0)|(a<=-1)|(b<=-1)) error('parameter(s) out of range'), end
sa=vpa(a); sb=vpa(b);
nu=(sb-sa)/(sa+sb+2);
mu=2^(sa+sb+1)*gamma(sa+1)*gamma(sb+1)/gamma(sa+sb+2);
if N==1, ab=[nu mu]; return, end 
N=N-1; n=1:N; nab=2*n+sa+sb;
A=[nu (sb^2-sa^2)*ones(1,N)./(nab.*(nab+2))];
n=2:N; nab=nab(n);
B1=4*(sa+1)*(sb+1)/((sa+sb+2)^2*(sa+sb+3));
B=4*(n+sa).*(n+sb).*n.*(n+sa+sb)./((nab.^2).*(nab+1).*(nab-1));
ab=[A' [mu; B1; B']];
