python compute the jacobian of numpy.array of function

Question

I want compute the numeric Jacobian of system of function contained into a numpy.array my variable are indexed in this way :

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

and then:

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

while my derivate is computed by :

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

the u[i] are the dependent variable while t is the independent one, I would find a way in order to compute the Jacobian using numerical differentiation (I want create a general case for a system of each type)

I find this one for matlab ... but I'm don't know it :(

function [J]=jacobian(func,x)
% computes the Jacobian of a function
n=length(x);
fx=feval(func,x);
eps=1.e-8; % could be made better
xperturb=x;
for i=1:n
xperturb(i)=xperturb(i)+eps;
J(:,i)=(feval(func,xperturb)-fx)/eps;
xperturb(i)=x(i);
end;