function [varargout]=lafunm(n,x)
  
% lafunm  Laguerre funcation of degree up to n defined by e^{-x/2}*lapolym(n,x)
   % y=lafunm(n,x) returns the Laguerre functions 
    % The degree should be a nonnegative integer, and the argument x is a vector 
    % [dy,y]=lafunm(n,x) also returns the values of 1st-order 
    %  derivatives of the Laguerre polynomials upto n stored in dy
    % Note: y (and likewise for dy) saves L_0(x), L_1(x), ...., L_n(x) by rows
    % i.e., L_k(x) is the (k+1)th row of the matrix y (or dy)
% Last modified on Decemeber 21, 2011        

dim=size(x); xx=x; if dim(1)>dim(2), xx=xx'; end; % xx is a row-vector     

if nargout==1,
   if n==0, varargout{1}=exp(-xx/2);  return; end;
   if n==1, varargout{1}=[exp(-xx/2);(1-xx).*exp(-xx/2)]; return; end;

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

if nargout==2,
   if n==0, varargout{2}=exp(-xx/2);  varargout{1}=-exp(-xx/2)/2; return; end;
   if n==1, varargout{2}=[exp(-xx/2);(1-xx).*exp(-xx/2)]; 
            varargout{1}=[-exp(-xx/2)/2; (xx-3).*exp(-xx/2)/2]; return; end;
   
   polylst=exp(-xx/2);  poly=(1-xx).*exp(-xx/2); pder=(xx-3).*exp(-xx/2)/2;
   y=[polylst;poly]; dy=[-exp(-xx/2)/2;pder];
   for k=1:n-1,
      polyn=((2*k+1-xx).*poly-k*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;
      y=[y;polyn]; dy=[dy; pdern]; 
   end;
      varargout{2}=y;  varargout{1}=dy; 
 end;
 
 return;