function [varargout]=hefunm(n,xx)
  
% hefunm  Hermite funcation of degree up to n defined by c_n e^{-x^2/2}*hepoly(n,x)
%        c_n=1/(pi^(1/4)(2^n n!)^(1/2)). 
    % y=hefunm(n,x) returns the Hermite functions 
    % The degree should be a nonnegative integer, and the argument x is a vector 
    % [dy,y]=hefunm(n,x) also returns the values of 1st-order 
    %  derivatives of the Hermite funcations up to 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 22, 2011       

   cst=1/sqrt(sqrt(pi));   
   dim=size(xx); x=xx; if dim(1)>dim(2), x=xx'; end; % x is a row-vector     
   
if nargout==1,    
   if n==0, varargout{1}=cst*exp(-x.^2/2);  return; end;
   if n==1, varargout{1}=[cst*exp(-x.^2/2);cst*sqrt(2)*x.*exp(-x.^2/2)]; return; end;

   polylst=cst*exp(-x.^2/2);	poly=cst*sqrt(2)*x.*exp(-x.^2/2);     
   y=[polylst;poly];
   for k=1:n-1,
	  polyn=sqrt(2/(k+1))*x.*poly-sqrt(k/(k+1))*polylst;  % see (7.73)
      polylst=poly; poly=polyn;	
      y=[y;poly];
   end;
   varargout{1}=y;
end;

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

  polylst=cst*exp(-x.^2/2);	poly=cst*sqrt(2)*x.*exp(-x.^2/2);     
  y=[polylst;poly]; dy=[cst*exp(-x.^2/2);cst*sqrt(2)*x.*exp(-x.^2/2)];
   for k=1:n-1,
      polyn=sqrt(2/(k+1))*x.*poly-sqrt(k/(k+1))*polylst;  % see (7.73)
      pdern=sqrt(2*(k+1))*poly-x.*polyn;    % see (7.75)      
      polylst=poly; poly=polyn;	      
      y=[y;poly]; dy=[dy;pdern];
   end;
      varargout{2}=y;  varargout{1}=dy;
 end;
 
 return