function [A,B]=jacobi2(fx_fil,t,x,u)
%JACOBI Function to calculate the Jacobian matrix numerically
%       [A,B]=jacobi2('fx_file',t,x,u)
%       PURPOSE
%       Function to calculate the Jacobian matrices A=df/dx and B=df/du
%       of the nonlinear function fx=f(t,x,u), i.e. the linearization of a possible non-linear
%       function fx=f(t,x,u) around values t x, and u.
%       ON INPUT
%       fx_file - m-file to define the function fx=f(t,x) with
%                 syntax: function fx=fx_file(t,x,u) 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
%       u  The input vector
%       ON OUTPUT
%       A       - The jacobian matrix A=df/dx of dimension (m,n)
%       B       - The Jacobian df/du of dimension (m=n,r)

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

uh=u;
B=zeros(m,nu);
for i=1:nu
   uh(i)=u(i)+h;
   fuh=feval(fx_fil,t,x,uh);
   uh(i)=uh(i)-h;
   B(:,i)=(fuh-fx)/h;
end
% END JACOBI2