function [pl,dpl]=hermite(n,x);

%       ==========================================================
%       Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
%                or Ln(x) or Hn(x), and their derivatives
%       Input :  KF --- Function code
%                       KF=1 for Chebyshev polynomial Tn(x)
%                       KF=2 for Chebyshev polynomial Un(x)
%                       KF=3 for Laguerre polynomial Ln(x)
%                       KF=4 for Hermite polynomial Hn(x)
%                n ---  Order of orthogonal polynomials
%                x ---  Argument of orthogonal polynomials
%       Output:  PL(n) --- Tn(x) or Un(x) or Ln(x) or Hn(x)
%                DPL(n)--- Tn'(x) or Un'(x) or Ln'(x) or Hn'(x)
%       =========================================================


kf = 4;

a=2.0d0;
b=0.0d0;
c=1.0d0;
y0=1.0d0;
y1=2.0d0.*x;
dy0=0.0d0;
dy1=2.0d0;
pl(0+1)=1.0d0;
pl(1+1)=2.0d0.*x;
dpl(0+1)=0.0d0;
dpl(1+1)=2.0d0;
if (kf == 1);
y1=x;
dy1=1.0d0;
pl(1+1)=x;
dpl(1+1)=1.0d0;
elseif (kf == 3);
y1=1.0d0-x;
dy1=-1.0d0;
pl(1+1)=1.0d0-x;
dpl(1+1)=-1.0d0;
end;
for  k=2:n;
if (kf == 3);
a=-1.0d0./k;
b=2.0d0+a;
c=1.0d0+a;
elseif (kf == 4);
c=2.0d0.*(k-1.0d0);
end;
yn=(a.*x+b).*y1-c.*y0;
dyn=a.*y1+(a.*x+b).*dy1-c.*dy0;
pl(k+1)=yn;
dpl(k+1)=dyn;
y0=y1;
y1=yn;
dy0=dy1;
dy1=dyn;
end;