function [nodes, weights] = GaussQuadrature(a,b,n)
%
% Returns the nodes and weights for the Gaussian quadrature on
% [-1,1] with weight function (1-x)^a * (1+x)^b. This is commonly called
% Gauss-Jacobi quadrature.
%

nmax = 50;

nodes  = zeros(n,1);
weights= zeros(n,1);
ab = a+b;

for i=1:n % Searching for node i
    
    % make an initial guess
    
    if (i==1) 
        an = a/n;
        bn = b/n;
        r1 = (1+a)*(2.78/(4+n*n)+0.768*an/n);
        r2 = 1+1.48*an+0.96*bn+0.452*an*an+0.83*an*bn;
        z  = 1-r1/r2;
    end;
    if (i==2) 
        r1 = (4.1+a)/((1+a)*(1+0.156*a));
        r2 = 1+0.06*(n-8)*(1+0.12*a)/n;
        r3 = 1+0.012*b*(1+0.25*abs(a))/n;
        z  = z - (1-z)*r1*r2*r3;
    end;
    if (i==3)
        r1 = (1.67+0.28*a)/(1+0.37*a);
        r2 = 1+0.22*(n-8)/n;
        r3 = 1+8*b/((6.28+b)*n*n);
        z  = z-(nodes(1)-z)*r1*r2*r3;  
    end;
    if ((i==n-1) && (i>3))
        r1 = (1+0.235*b)/(0.766+0.119*b);
        r2 = 1/(1+0.639*(n-4)/(1+0.71*(n-4)));
        r3 = 1/(1+20*a/((7.5+a)*n*n));
        z  = z+(z-nodes(n-3))*r1*r2*r3;
    end;
    if ((i==n) && (i>3))
        r1 = (1.0+0.37*b)/(1.67+0.28*b);
        r2 = 1/(1+0.22*(n-8)/n);
        r3 = 1/(1+8*a/((6.28+a)*n*n));
        z  = z+(z-nodes(n-2))*r1*r2*r3;
    end;
    
    if ((i>3) && (i<n-1))
        z = 3*nodes(i-1)-3*nodes(i-2)+nodes(i-3);
    end;
        
    its   = 1;
    found = 0;
    while ( (its<=nmax) && (found==0))
       temp = 2+ab;
       p1 = (a-b+temp*z)/2;
       p2 = 1;
       for j=2:n % evaulate the polynomial
           p3 = p2;
           p2 = p1;
           temp  = 2*j+ab;
           temp1 = 2*j*(j+ab)*(temp-2);
           temp2 = (temp-1)*(a*a-b*b+temp*(temp-2)*z);
           temp3 = 2*(j-1+a)*(j-1+b)*temp;
           p1    = (temp2*p2-temp3*p3)/temp1;
       end;
       
       % and its derivative
       pp = (n*(a-b-temp*z)*p1+2*(n+a)*(n+b)*p2)/(temp*(1-z*z));
       
       % Newton step
       z1 = z;
       z  = z1-p1/pp;
       
       its = its+1;
       if (abs(z1-z)<1e-14)
           found = 1;
       end;
    end;
    
    if (found==0) 
        disp('Gaussian quadrature points not found !');
    else
        nodes(i) = z;
        weights(i) = exp(gammaln(a+n)+gammaln(b+n)-gammaln(n+1)-gammaln(n+ab+1))...
                     *temp*2^(ab)/(pp*p2);
    end
    
end % GaussQuadrature