function [varargout]=glagsrd(n,alp)

% x=glagsrd(n,alp) returns n generalized Laguerre-Gauss-Radau points with x(1)=0.
%  [x,w]= glagsrd(n,alp)  returns n  Generalized Laguerre-Gauss-Radau 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

% Compute the zeros of L_{n}^{(alp)}'(x) i.e., L_{n-1}^{(alp+1)}(x)
J=diag(2*[0:n-2]+alp+2)+diag(-sqrt([1:n-2].*([1:n-2]+alp+1)),1)+diag(-sqrt([1:n-2].*([1:n-2]+alp+1)),-1);  % (7.38) of the book
r = sort(eig(sparse(J)));               % Compute eigenvalues

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

% Compute the weights (7.28) of the book 
omega0=(alp+1)*(gamma(alp+1))^2*exp(gammaln(n)-gammaln(n+alp+1));  % omega_0
y=glapoly(n-1,alp,r);                 
gm=gammaln(n+alp)-log(n+alp)-gammaln(n); gm=exp(gm);
varargout{2}=[omega0; gm./(y.^2)];                                  % omega_j

if nargout==2, return; end;

% Compute the weights of  generalized Laguerre function  approach  (7.34)  of the book
z=glafun(n-1,alp,r);
varargout{3}=[omega0; gm./(z.^2)];  

return