function fdg = fd_grad(fname,u)
% Usage: fdg = fd_grad(fname,u)
% Purpose: computes the finite difference gradient for function fname at u
% Input:
%   fname: name of the function to be evaluated
%   u: point where the finite difference gradient is to be computed
% Output:
%   fdg: finite difference approximation of the gradient of fname at u
% Author: TK 3.11.2002 (also late in the dark evening...)
delta = sqrt(eps); %FD stepsize: square root of machine epsilon
n = length(u);
fdg = zeros(size(u));
for i=1:n
  v = u;
  v(i) = v(i) + delta;
  fpd = feval(fname,v);
  v = u;
  v(i) = v(i) - delta;
  fmd = feval(fname,v);
  fdg(i) = (fpd - fmd)/(2*delta);
end%for
