function [x,w]=gauss_jacobi(N,alpha,beta)

% GAUSS-JACOBI RULE ON [-1,1] (OR (-1,1) DEPENDING ON alpha, beta).
% N: FOR GAUSSIAN RULES IT IS THE NUMBER OF QUAD. POINTS. 
% alpha, beta ARE THE GAUSS-JACOBI EXPONENTS.
% x, w ARE COLUMN VECTORS OF NODES AND WEIGHTS. 
%      THE LENGTH OF x AND w IS "N" IF gl=0, "N+2" IF "gl=1".

if nargin < 2
    alpha=0; beta=0;
end

ab=r_jacobi(N,alpha,beta);
xw=gauss(N,ab);
    
x=xw(:,1);

w=xw(:,2);

%--------------------------------------------------------------------------
% ADDITIONAL FUNCTIONS BY D.LAURIE AND W.GAUTSCHI.
%--------------------------------------------------------------------------

function ab=r_jacobi(N,a,b)

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];


function xw=gauss(N,ab)
N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
  J(n,n-1)=sqrt(ab(n,2));
  J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];
