def METHOD_NAME(x):
    v0 = f(x)
    N = v0.shape[0]
    n_v = v0.shape[1]
    assert x.shape[0] == N
    n_x = x.shape[1]
    dv = zeros((N, n_v, n_x))
    for i in range(n_x):
        xi = x.copy()
        xi[:, i] += epsilon
        vi = f(xi)
        dv[:, :, i] = (vi - v0) / epsilon
    return [v0, dv]