function [varargout]=glags(n,alp)

% x=glags(n,alp) returns n generalized Laguerre-Gauss points 
%  [x,w]= glags(n,alp)  returns n  Generalized Laguerre-Gauss points and weights
%  [x,w,wf]=glagsrd(n,alp) also returns the weights (in wf) associated with generalized Laguerre
%  function approach (see (7.34) of the book). 
%  Eigenmethod is used for computing nodes.
% Use: glapoly( ); glafun();
%  Last modified on Decemeber 21, 2011


J=diag(2*[0:n-1]+alp+1)+diag(-sqrt([1:n-1].*([1:n-1]+alp)),1)+diag(-sqrt([1:n-1].*([1:n-1]+alp)),-1);  % (7.38) of the book
r = sort(eig(sparse(J)));               % Compute eigenvalues

varargout{1}=r;
if nargout==1, return; end;

% Compute the weights (7.27) of the book 
y=glapoly(n-1,alp,r);
gm=gammaln(n+alp)-log(n+alp)-gammaln(n+1);gm=exp(gm);
varargout{2}=gm*r./(y.^2);  

if nargout==2, return; end;

z=glafun(n-1,alp,r);
varargout{3}=gm*r./(z.^2);  % output weights related to Laguerre function approach.

return