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