function A=jacobi(fx_fil,t,x)
%JACOBI Function to calculate the Jacobian matrix numerically
%       A=jacobi('fx_file',t,x)
%       PURPOSE
%       Function to calculate the Jacobian matrix A=df/dx of the nonlinear
%       function fx=f(t,x), i.e. the linearization of a possible non-linear
%       function fx=f(t,x) around values t and x.
%       ON INPUT
%       fx_file - m-file to define the function fx=f(t,x) with
%                 syntax: function fx=fx_file(t,x) where fx is, 
%                 say  m dimensional
%       t       - time instant.
%       x       - column vector with same dimension as in the function
%           fx=f(t,x), say of dimension n
%       ON OUTPUT
%       A       - The jacobian matrix A=df/dx of dimension (m,n)

h=1e-5;      % Fixed step length in dot(f) = (f(x+h)-f(x))/h
n=length(x);
xh=x;
fx=feval(fx_fil,t,x);
m=length(fx);
A=zeros(m,n);
for i=1:n
   xh(i)=x(i)+h;
   fxh=feval(fx_fil,t,xh);
   xh(i)=xh(i)-h;
   A(:,i)=(fxh-fx)/h;
end
% END JACOBI
