function [x,w] = legendre(a,b,n);
%
% [x,w] = legendre(a,b,n)
%
% n-point Gauss quadrature rule on [a,b] with weight function
% w(t) = 1/(b-a), a <= t <= b.
%
% INPUT
%
% a,b  : left and right interval end point
% n    : number of nodes, n > 0
%
% OUTPUT
%
% x   : n x 1 vector containing cubature nodes
% w   : n x 1 vector containing cubature weights
%
% \int_a^b f(t) * w(t) dt = \sum_{n=1}^n f(x(n,:))*w(n)
%
% Copyright (c) 2007 E. Ullmann
%

alpha = zeros(n,1);
m = [2:n]';
beta = (m-1) ./ sqrt((2*m-3).*(2*m-1));

% Rescale recurrence coefficients to interval [a,b]
        
beta = (b-a)/2 * beta;
alpha = a + (b-a)/2 * (1 + alpha);
        
% Set up Jacobi matrix and compute eigendecomposition
%
J = spdiags(beta,-1,n,n);
J = J + J' + spdiags(alpha,0,n,n);
    
[C,mu] = eig(full(J));
x = diag(mu);
w = (C(1,:)').^2;

return