Python Compute Jacobian numerically Hi I have a class with the purpose to solve system of Differential equations, this class contains a class named Rhs (right hand side) that contains all the feature of the differential problem: init_time , final_time , Initial value , function(that is a numpy.array containing lambda function)
when I have a system of differential equation the functions are managed in this way :
eq1 = lambda t,u : a*(u[0]-u[0]*u[1]);
eq2 = lambda t,u : -c*(u[1]-u[0]*u[1]);

func1 = np.array([eq1,eq2])

y0      = np.array([2.,1.])
system1 = rhs.Rhs(func1, 0,10,y0,500 )

using the parameter above the class Rhs is completed ! it have a method to return the value of each function, and a method for compute the derivate of each of function:
def f(self,ti,ui):
       return  np.array([function(ti,ui) for function in self.func])     


def Df(self,ti,ui):
      eps = 10e-12
      return (self.f(ti,ui+eps) - self.f(ti,ui-eps) )/(2*eps)

Now I have a big problem, in order to create a class to solve an implicit method I have to compute the Jacobian of the function ! but I have no idea how to do this ! 
EDIT no I need to define the jacobian matrix yes I wrote the method for derivative but I really have not idea how to define J[i][j] @saulspatz may you help me pls ? I don't need to compute it's determinant just to write the terms $J_{i,j}$
I wrote this but doesn't works :
   def J(self,ti,ui):   
     self.t = ti 
     self.u = ui
     Jac = np.zeros((len(ui),len(ui)))
     n = len(ui)
     eps = 1e-12
     for i in range(n):
         for j in range(n):
             Jac[i,j] = self.Df(ti,ui) 

     return Jac

EDIT thank you very much @caverac how can I generalize for any system ? 
is that straightforward ? 
 A: If you have a problem of the form
\begin{eqnarray}
\frac{{\rm d}x_1}{{\rm d}t} &=& f_1(x_1, x_2, \cdots) \\
\frac{{\rm d}x_2}{{\rm d}t} &=& f_2(x_1, x_2, \cdots) \\
&\vdots&
\end{eqnarray}
then the Jacobian you are looking for is
$$
J = \pmatrix{
\partial f_1/\partial x_1 & \partial f_1/\partial x_2 & \cdots \\
\partial f_2/\partial x_1 & \partial f_2/\partial x_2 & \cdots \\ 
\vdots & \vdots & \ddots
}
$$
or
$$
J_{ij} = \frac{\partial f_i}{\partial x_j}
$$
in your case
$$
f_1(x_1, x_2) = a(x_1 - x_1 x_2) ~~~\mbox{and}~~ f_2(x_1, x_2) = -c (x_2 - x_1 x_2)
$$
So the Jacobian takes the form
$$
J = \pmatrix{
a(1 - x_2) & -ax_1 \\
cx_2 & -c(1 - x_1)
}
$$
def f(t, x, **params):

    a = params['a']
    c = params['c']

    f1 = a * (x[0] - x[0] * x[1])
    f2 = -c * (x[1] - x[0] * x[1])   

    return np.array([f1, f2], dtype = np.float)

def df(t, x, **params):

    eps = 1e-10
    J = np.zeros([len(x), len(x)], dtype = np.float)

    for i in range(len(x)):
        x1 = x.copy()
        x2 = x.copy()

        x1[i] += eps
        x2[i] -= eps

        f1 = f(t, x1, **params)
        f2 = f(t, x2, **params)

        J[ : , i] = (f1 - f2) / (2 * eps)

    return J

t = 0
x = np.array([1, 2], dtype = np.float)
print df(t, x, a = 1, c = 1)

EDIT Wrapped methods into a class
class ODEModel(object):

    def __init__(self, eps = 1e-10):

        self.eps = eps
        self.farray = []

    def add_function(self, f):

        self.farray.append(f)

    def f(self, t, x):

        return np.array([fi(t, x) for fi in self.farray], dtype = np.float)

    def df(self, t, x):

        J = np.zeros([len(x), len(x)], dtype = np.float)

        for i in range(len(x)):
            x1 = x.copy()
            x2 = x.copy()

            x1[i] += self.eps
            x2[i] -= self.eps

            f1 = self.f(t, x1)
            f2 = self.f(t, x2)

            J[ : , i] = (f1 - f2) / (2 * self.eps)

        return J


F = ODEModel(eps = 1e-12)

eq1 = lambda t,u : u[1]
eq2 = lambda t,u : u[2]
eq3 = lambda t,u : u[3]
eq4 = lambda t,u : -8 * u[0] + np.sin(t) * u[1] - 3 * u[3] + t**2

F.add_function(eq1)
F.add_function(eq2)
F.add_function(eq3)
F.add_function(eq4)


t = 0
x = np.array([1, 2, 3, 4], dtype = np.float)
print F.df(t, x)

