function f = hermfunc(m, y)

if m == 0,
   f = exp(-0.5 * y.^2) / (pi .^ 0.25);
elseif m == 1,
   f = (2/sqrt(2)) * y .* exp(-0.5 * y.^2) / (pi .^ 0.25);
else
   f = sqrt(2/m) * y .* hermfunc(m-1, y) ...
         - sqrt(1 - (1/m)) * hermfunc(m-2, y);
end
