function ab=r_jacobi(N,a,b)
% This fiunction returns the recurrence coefficients for monic Jacobi polynomials.
%
%    ab=r_jacobi(N,a,b) generates the first N recurrence 
%    coefficients for monic Jacobi polynomials with parameters 
%    a and b. These are orthogonal on [-1,1] relative to the
%    weight function w(t)=(1-t)^a(1+t)^b. The n alpha-coefficients
%    are stored in the first column, the n beta-coefficients in
%    the second column, of the nx2 array ab.
% 
% References:
% Walter Gautschi website
%  
% INPUTS
% N - number of nodes
% a,b - esponenti della funzione w(x)=((1-x)^a)*((1+x)^b);
% 
% OUTPUTS
% ab - the n alpha-coefficients are stored in the first column,
% the n beta-coefficients in the second column, of the nx2 array
% abcoefficienti di ricorrenza dei polinomi ortogonali

nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1
 ab=[nu mu]; return
end

N=N-1;
n=1:N;
nab=2*n+a+b;
nuadd=(b^2-a^2)*ones(1,N)./(nab.*(nab+2));
A=[nu nuadd];
n=2:N;
nab=nab(n);
B1=4*(a+1)*(b+1)/((a+b+2)^2*(a+b+3));
B=4*(n+a).*(n+b).*n.*(n+a+b)./((nab.^2).*(nab+1).*(nab-1));
abadd=[mu; B1; B'];
ab=[A' abadd];