function retval = numhess(fname,x)
% usage:  y = numhess(fname,x)
% description: returns numerical function Hessian
% using forward differences

global parms;

% no bulletproofing
n = length(x);
h = parms.hstep;
h2 = h*h;
hvec = h*eye(n);
fx = feval(fname,x);
fxdelta = zeros(n,1);
retval = zeros(n,n);
for jj = 1:n
  fxdelta(jj) = feval(fname,x+hvec(:,jj));
end
for jj = 1:n
  for kk = jj:n
    retval(jj,kk) = (feval(fname,x+hvec(:,jj)+hvec(:,kk)) + fx - ...
           fxdelta(jj) - fxdelta(kk))/h2;
  end
end
retval = retval + triu(retval,1)';