def approx_jacobian(x,func,epsilon,*args):
x0 = asfarray(x)<EOL>f0 = atleast_1d(func(*((x0,)+args)))<EOL>jac = zeros([len(x0),len(f0)])<EOL>dx = zeros(len(x0))<EOL>for i in range(len(x0)):<EOL><INDENT>dx[i] = epsilon<EOL>jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon<EOL>dx[i] = <NUM_LIT:0.0><EOL><DEDENT>return jac.transpose()<EOL>
Approximate the Jacobian matrix of a callable function.

Parameters
----------
x : array_like
    The state vector at which to compute the Jacobian matrix.
func : callable f(x,*args)
    The vector-valued function.
epsilon : float
    The perturbation used to determine the partial derivatives.
args : sequence
    Additional arguments passed to func.

Returns
-------
An array of dimensions ``(lenf, lenx)`` where ``lenf`` is the length
of the outputs of `func`, and ``lenx`` is the number of elements in
`x`.

Notes
-----
The approximation is done using forward differences.