function J=jacobian(f,x,t)
  f0=f(x,t);
  J=zeros(length(x),length(x))
  for i=1:length(x)
      dx=zeros(x);
      dx(i)=abs(x(i))*1e-6+1e-9;
      xp=x+dx;
      J(:,i)=(f(xp,t)-f0)/dx(i);
  end
endfunction
