def jacobian(ctx, f, x):
    x = ctx.matrix(x)
    h = ctx.sqrt(ctx.eps)
    fx = ctx.matrix(f(*x))
    m = len(fx)
    n = len(x)
    J = ctx.matrix(m, n)
    for j in xrange(n):
        xj = x.copy()
        xj[j] += h
        Jj = (ctx.matrix(f(*xj)) - fx) / h
        for i in xrange(m):
            J[i,j] = Jj[i]
    return J

# TODO: test with user-specified jacobian matrix, support force_type
Calculate the Jacobian matrix of a function at the point x0.
dataset/ETHPy150Open fredrik-johansson/mpmath/mpmath/calculus/optimization.py/jacobian