function J=NumericalJacobian(f,x);
% NUMERICALJACOBIAN computes a Jacobian numerically
%   J=NumericalJacobian(f,x); computes numerically a Jacobian matrix
%   for the function f at x
    
for i=1:length(x)
  xd=x; h=sqrt(eps*(1+abs(xd(i)))); xd(i)=xd(i)+h; 
  J(:,i)=(f(xd)-f(x))/h;
end;