function [x,w] = GaussLaguerre(n,alpha)
%      [x,w] = GaussLaguerre(n,alpha)
%      find Gauss-Laguerre nodes x(1),...x(n) and weights w(1),...,w(n)
%      need alpha>-1
%
%      integral of f(x)*x^alpha*exp(-x) from 0 to infinity
%      is approximated by w(1)*f(x(1)) + ... + w(n)*f(x(n))
%
%      integral of f(x)*x^alpha*exp(-b*x) from 0 to infinity    (with b>0)
%      is approximated by b^(-1-alpha)*(w(1)*f(t(1)) + ... + w(n)*f(t(n)))
%      where t=x/b
%
% Example: find integral of cos(3*x)*x^(-1/3)*exp(-2*x) from 0 to infinity
%   alpha = -1/3;
%   [x,w] = GaussLaguerre(20,alpha);
%   b=2; t=x/b;
%   Q = b^(-1-alpha)*sum(w.*cos(3*t))
 
% from www.mathworks.com/matlabcentral/fileexchange/26737-legendre-laguerre-and-hermite-gauss-quadrature
%  Geert Van Damme
% geert@vandamme-iliano.be
% February 21, 2010    

% Building the companion matrix CM
    % CM is such that det(xI-CM)=L_n(x), with L_n the Laguerree polynomial
    % under consideration. Moreover, CM will be constructed in such a way
    % that it is symmetrical.
i   = 1:n;
a   = (2*i-1) + alpha;
b   = sqrt( i(1:n-1) .* ((1:n-1) + alpha) );
CM  = diag(a) + diag(b,1) + diag(b,-1);

% Determining the abscissas (x) and weights (w)
    % - since det(xI-CM)=L_n(x), the abscissas are the roots of the
    %   characteristic polynomial, i.d. the eigenvalues of CM;
    % - the weights can be derived from the corresponding eigenvectors.
[V L]   = eig(CM);
[x ind] = sort(diag(L));
V       = V(:,ind)';
w       = gamma(alpha+1) .* V(:,1).^2;
