def METHOD_NAME(x0, fn, epsilon):
    """Computes the Jacobian using finite differences
    x0:      The positions at which to compute J.
    fn:      A function of the form fn(x) which returns a single numpy array.
    epsilon: A scalar or an array that can be broadcasted to the same
             shape as x0.
    """
    dtype = x0.dtype
    y0 = fn(x0)
    h = np.zeros_like(x0)
    J = np.zeros((x0.size, y0.size), dtype=dtype)
    epsilon_arr = np.broadcast_to(epsilon, x0.shape)
    for i in range(x0.size):
        eps = epsilon_arr.flat[i]
        h.flat[i] = eps
        J[i, :] = ((fn(x0 + h) - y0) / eps).flat
        h.flat[i] = 0
    return J