function retval = numgrad(fname,x)
% usage:  y = numgrad(fname,x)
% description: returns numerical function gradient
% using central differences

global parms;

% no bulletproofing
n = length(x);
h = parms.hstep;
twoh = 2*h;
retval = zeros(size(x));
for jj = 1:n
  xup = x;
  xup(jj) = xup(jj)+h;
  xdown = x;
  xdown(jj) = xdown(jj)-h;
  retval(jj) = (feval(fname,xup)-feval(fname,xdown))/twoh;
end