% hermite functions
function y = hermiteS(N, x, scale, shift) % works well

if (nargin < 4)
    shift = 0;
end
if (nargin < 3)
    scale = 1;
end
if (N ~= round(N))
    error('??? N has to be integer valued.');
end
if (~isscalar(N))
    si = size(x);
    if (si(2) == 1)
        x = transpose(x);
    end
    y = zeros(length(N), 1) * x;
end

y0 = 0 * x;
x = scale * (x - shift);
y1 = (scale^2 / pi)^(0.25) * exp(- x.^2 / 2);

k1 = 1;
if (N(1) == 0)
    y(1,:) = y1;
    k1 = k1+1;
end

for k=1:max(N)
    y2 = sqrt(2.0 / k) * x .* y1 - sqrt(1 - 1.0 / k) * y0;
    y0 = y1;
    y1 = y2;

  if (k == N(k1))
      y(k1,:) = y1;
      k1 = k1 + 1;
  end
end

return