function [varargout]=lapolym(n,x)

% lapolym  Laguerre polynomials of degree up to n, i.e., 0,1, ...,n
    % y=lapolym(n,x) returns the Laguerre polynomials 
    % The degree should be a nonnegative integer, and the argument x is a vector 
    % [dy,y]=lapolym(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}=ones(size(xx));  return; end;
   if n==1, varargout{1}=[ones(size(xx));1-xx]; return; end;

   polylst=ones(size(xx));	poly=1-xx;  % L_0=1; L_1=1-x;    
   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}=ones(size(xx));  varargout{1}=zeros(size(xx)); return; end;
   if n==1, varargout{2}=[ones(size(xx));1-xx]; 
            varargout{1}=[zeros(size(xx)); -ones(size(xx))]; return; end;

   polylst=ones(size(xx)); poly=1-xx; pder=-ones(size(xx)); % L_0=1, L_1=1-x, L_1'=-1; 
   y=[polylst;poly]; dy=[zeros(size(xx));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-poly ;	                	   % L_{k+1}'=L_k'-L_k; 
	  polylst=poly;  poly=polyn;  pder=pdern; 
      y=[y;polyn]; dy=[dy; pdern]; 
   end;
      varargout{2}=y;  varargout{1}=dy;
 end;
 
 return;