function J = appr_dfdx(f,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% function J = appr_dfdx(f,x) %
% %
% Numerical approximation of f's jacobian %
% at x. f is a string with the name of the %
% function file in which f is defined %
% %
% Uses the simple approximation %
% f'(x) = (f(x+h)-f(x))/h %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = length(x);
fx = feval(f,x);
m = length(fx);
h = 1E-4;
for i = 1:n
xph = x;
xph(i) = x(i)+h;
fxph = feval(f,xph);
J(:,i) = (fxph-fx)/h;
end
err = 1;
tol = 1E-2;
while (err>tol)
h_old = h;
J_old = J;
h = h/2;
for i = 1:n
xph = x;
xph(i) = x(i)+h;
fxph = feval(f,xph);
J(:,i) = (fxph-fx)/h;
end
err = norm(J-J_old,'inf');
errorInDfDx=err; % (trace output)
end