Func Matrix numerical_jacobi(f, t, x, xeps)
     // ヤコビ行列を線形近似で求める
     Matrix f();
     Real t;
     Matrix x, xeps;
{
     Integer dim, k;
     Matrix jacobian;
     Matrix epsm;

     dim= length(x);
     epsm = vec2diag(xeps);

/*
  jacobian = 

 +-                                                -+
 | jac11=df1/dx1  jac12=df1/dx2  jac13=df1/dx3  ... |
 | jac21=df2/dx1  jac22=df2/dx2  jac23=df2/dx3  ... |
 | jac31=df3/dx1  jac32=df3/dx2  jac33=df3/dx3  ... |
 | ............                                     |
 +-                                                -+
*/

     jacobian = Z(dim);
     for (k=1; k<=dim; k++){
	  jacobian(:,k)= Array( f(t,x+epsm(:,k))-f(t,x) )/Array(xeps(k)); // 各列
     }
     return jacobian;
}