function [varargout]=glafun(n,alp,x)
  
% glafun  generalized Laguerre funcation of degree n defined by e^{-x/2}*glapoly(n,alp,x)
    % y=glafun(n,alp,x) returns the generalized Laguerre function of degree n at x
    % The degree should be a nonnegative integer 
    % The arguments x >0, alp>-1; 
    % [dy,y]=glafun(n,alp,x) also returns the first-order 
    %  derivative of the generalized Laguerre function in dy
% Last modified on Decemeber 21, 2011        
     
if nargout==1,
   if n==0, varargout{1}=exp(-x/2);  return; end;
   if n==1, varargout{1}=(alp+1-x).*exp(-x/2); return; end;

   polylst=exp(-x/2);	poly=(alp+1-x).*exp(-x/2);  % L_0=e^{-x/2}; L_1=(1-x)*e^{-x/2};    
   for k=1:n-1,
	  polyn=((2*k+alp+1-x).*poly-(k+alp)*polylst)/(k+1);  % L_{k+1}=((2k+1-x)L_k-kL_{k-1})/(k+1);
      polylst=poly; poly=polyn;	
   end;
   varargout{1}=poly;    
end;

if nargout==2,
   if n==0, varargout{2}=exp(-x/2);  varargout{1}=-exp(-x/2)/2; return; end;
   if n==1, varargout{2}=(alp+1-x).*exp(-x/2); varargout{1}=-(alp+3-x).*exp(-x/2)/2; return; end;
   
   polylst=exp(-x/2);  poly=(alp+1-x).*exp(-x/2); pder=-(alp+3-x).*exp(-x/2)/2;
   for k=1:n-1,
      polyn=((2*k+alp+1-x).*poly-(k+alp)*polylst)/(k+1); % L_{k+1}=((2k+1-x)L_k-kL_{k-1})/(k+1);
	  pdern=pder-(polyn+poly)/2;	       	   % L_{k+1}'=L_k'-(L_{k+1}+L_k)/2;  by (7.26b) of the book 
	  polylst=poly; poly=polyn; pder=pdern;
   end;
      varargout{2}=poly;  varargout{1}=pder;
 end;
 
 return;