function A=jacobi(f,x)
% Computes the derivative (Jacobi matrix) A=Df(x) of the function f at x.
% 
%   Syntax: 
%           A = jacobi(f,x)
%   Arguments: 
%           f - string containing the name of a function file  
%               which returns an mx1 matrix y=f(x).  
%           x - an nx1 matrix
%   Returns: 
%           A - an mxn matrix, the Jacobian of f 
%   Description:
%           The program jacobi uses a central difference quotient
%           to compute an approximation of the Jacobi matrix A=Df(x)
%           The file f.m must return an mx1 matrix y=f(x).  
%           The stepsize is h=1e-5.  See Body and Soul, Sect. 23.13.
%           The optimal stepsize for the central difference quotient is 
%           h = (df/Kf)^(1/3), where df=1e-16 is the error in f(x) and 
%           Kf is the constant in Taylor's formula, Sect. 24.8, 
%           f(x)=f(a)+f'(a)(x-a)+0.5 f''(a)(x-a)^2 + Ef, |Ef| < Kf|x-a|^3.
%           The error is then 2(df)^(2/3)(Kf)^(1/3), approximately 1e-10.
%   Examples: 
%          z = jacobi('sin', 0) returns 1
%          z = jacobi('funk1', 3) returns 6
%          A = jacobi('funk2', [1;2;3]) returns [2 4 6]
%   See also:
%          bisec.m, fixpoint.m, regula_falsi.m, newton2.m
%

%----------------------------------------------------------------------
%   Date created: 1999-12-05
%
%   Author:  Stig Larsson
%------------------------------+---------------------------------------+

h = 1e-5;   % stepsize in the central difference approximation

% test if f is a column vector
y=feval(f,x);
if size(y,2)~=1 
disp('jacobi: f is not a column vector'); 
return;
end

for i = 1:size(x,1)

   e = zeros(size(x));
   e(i) = 1;                 % e is the i:th basis vector

   % compute i:th column of A 
   A(:,i) = ( feval(f,x+h*e) - feval(f,x-h*e) )/(2*h);

end


