def METHOD_NAME(self):
    fun = lambda x: [-josephy(x), -Djosephy(x)]
    x0 = np.array([1.25, 0.01, 0.01, 0.50])
    lb = np.array([0.00, 0.00, 0.00, 0.00])
    # ub=np.array( [inf, inf, inf, inf] )
    ub = np.array([1e20, 1e20, 1e20, 1e20])
    # resp = ncpsolve(fun,  lb, ub, x0, tol=1e-15)
    sol = np.array([1.22474487e00, 0.00000000e00, 3.60543164e-17, 5.00000000e-01])
    # assert_almost_equal(sol,  resp)
    N = 10
    d = len(x0)
    s_x0 = np.row_stack([x0] * N)
    s_lb = np.row_stack([lb] * N)
    s_ub = np.row_stack([ub] * N)
    def serial_fun(xvec):
        resp = np.zeros((N, d))
        dresp = np.zeros((N, d, d))
        for n in range(N):
            [v, dv] = fun(xvec[n, :])
            resp[n, :] = v
            dresp[n, :, :] = dv
        return [resp, dresp]
    res = serial_fun(s_x0)[0]
    serial_sol = ncpsolve(
        serial_fun, s_lb, s_ub, s_x0, jactype="serial", verbose=True
    )
    print(serial_sol)