function J=jacobian(f,x)
  //calculates the jacobian dx/df evaluated at x
  n=length(x);
  J=zeros(n,n);
  fx=f(x);
  for i=1:n
    x1=x;
    dx=0.00001*x1(i)+1e-10;
    x1(i)=x1(i)+dx;
    J(:,i)=(f(x1)-fx)/dx;
  end
endfunction  
  
