In a general Gaussian quadrature rule, a definite integral
of f(x) is first approximated over the interval [-1, 1] by a
polynomial approximatedfunction g(x) and a known
weighting  function W(x).

$\int_{-1}^{1}f(x)dx=\int_{-1}^{1}W(x)g(x)dx$

Those are then approximated by a sum of function values
at  specified points $x_{i}$ multiplied by some weights w_i:

$\int_{-1}^{1}W(x)g(x)dx\approx \sum_{i=1}^{n}w_{i}g(x_{i})$

In the case of Gauss-Legendre quadrature, the weighting  
function W(x) = 1, so we can approximate an integral of
 f(x)  with:

$\int_{-1}^{1}f(x)dx\approx \sum_{i=1}^{n}w_{i}f(x_{i})$

The n evaluation points x_i for a  n-point rule, also
called "nodes", are roots of n-th order Legendre  
Polynomials P_n(x). Legendre polynomials are
defined  by the following recursive rule:

P₀(x)  = 1

P₁(x)  = x
nP_n(x) = (2n − 1)xP_{n −  1}(x) − (n − 1)P_{n − 2}(x)

There is also a recursive equation for their
derivative:

${P_{n}}'(x)=\frac{n}{x^{2}-1}\left [ xP_{n}(x)-P_{n-1}(x) \right ]$

The roots of those polynomials are in general not
analytically solvable, so they have to be
approximated numerically, for  example by
Newton-Raphson iteration:

$x_{n+1}=x_{n}-\frac{f(x_{n})}{{f}'(x_{n})}$

The first guess x₀ for the i-th  root of a n-order
polynomial P_n can be given by

$x_{0}=cos(\pi \frac{i-0.25}{n+0.5})$

After we get the nodes x_i, we compute  the
appropriate weights by:

  $w_{i}=\frac{2}{(1-x_{i}^{2})[{P}'_{n}(x_{i})]^{2}}$

After we have the nodes and the weights for a
n-point  quadrature rule, we can approximate an
integral over any interval [a,b]  by

$\int_{a}^{b}f(x)dx\approx \frac{b-a}{2}\sum_{i=1}^{n}w_{i}f(\frac{b-a}{2}x_{i}+\frac{b+a}{2})$

     function [w,x]=lgwt(N)            

   x=cos(pi*((1:N)'-0.25)/(N+0.5));    % Initial guess
   L=zeros(N,N+1);                     % Legendre-Gauss Matrix

   x0=2;
   while max(abs(x-x0))>eps
        L(:,1)=1;
        L(:,2)=x;
   
         for k=2:N
               L(:,k+1)=((2*k-1)*x.*L(:,k)-(k-1)*L(:,k-1))/k; % Legendre polynormial Pn defined by recursive rule
         end
   
        Lp=N*(x.*L(:,N+1)-L(:,N))./(x.^2-1);       % Derivative Pn'
   
        x0=x;
        x=x0-L(:,N+1)./Lp;          % Newton-Raphson method for the roots
   
   end

           w=2./((1-x.^2).*Lp.^2);         % Compute the weights