function [x, w] = gauss(N)
%   GAUSS(N) returns N Legendre points X in (-1,1).
%
%   [X, W] = GAUSS(N) returns also a row vector W of weights for
%   Gauss-Legendre quadrature.  If f is defined as an @ function
%   the integral of f on [-1,1] can be computed as sum(w.*f(x)).

% Nick Hale - Apr 2015

% Note: for simplicity we use the Chebyshev nodes, but one can do better.
x = -cos(((.5:N-.5).')*pi/N);   % Initial guess for roots of P_N(x).

tol = eps;                      % Tolerance (can be changed).
dx = inf;                       % Initialise convergence test.

% Loop until convergence:
while ( norm(dx, inf) > .1*sqrt(tol) )
    
    % Recurrence relation for Legendre polynomials:
    Pm2 = 1; Pm1 = x;                       % Initialise.
    for n = 1:(N-1)
        P = ( (2*n+1)*Pm1.*x - n*Pm2 ) / (n+1);
        Pm2 = Pm1; Pm1 = P;                 % Store.
    end
    dPdx = -N * (x.*P - Pm2) ./ (1 - x.^2); % See NIST (18.9.17).
    
    % Newton step:
    dx = -P ./ dPdx;
    % Newton update:
    x = x + dx;

end

% Quadrature weights:
w = 2 ./ ( (1-x.^2) .* dPdx.^2 );

end
