function ret = jacob(fun, x)
% function ret = jacob(fun, x)
%
% numerical differention of a function using symmetric
% differences method
%
% inputs:
%  fun - function whose derivative is to be computed;
% 	it is assumed that function accepts and returns column vector
%  x - point at which the derivative is to be computed (column vector)
%
% outputs:
%  ret - numerical derivative obtained
%
% (c) Grzegorz Klima 2004


m = size(x, 1);
n = size(feval(fun, x), 1);

h = eps ^ (1 / 3);

jac = zeros(n, m);

for j = (1:m)
  x1 = x;
  x2 = x;
  x1(j) = x(j) + h;
  x2(j) = x(j) - h;
  f1 = feval(fun, x1);
  f2 = feval(fun, x2);
  jac(:, j) = (f1 - f2) / (x1(j) - x2(j));
end

ret = jac;