import numpy as np
from scipy import optimize

J = 100.
N = 5000.
U = 0.05

Lambda = N*U/(2*J)
print Lambda, np.sqrt(1-1/(Lambda**2))


def func(x):
    return [-2 * J * np.sqrt(1 - x[0] ** 2) * np.sin(x[1]),
            N * U * x[0] + 2 * J * x[0] * np.cos(x[1]) / np.sqrt(1 - x[0] ** 2)]


def jac(x):
    return np.array([[2 * J * x[0] * np.sin(x[1]) / np.sqrt(1 - x[0] ** 2),
                      -2 * J * np.sqrt(1 - x[0] ** 2) * np.cos(x[1])],
                     [N * U + 2 * J * np.cos(x[1]) / (np.sqrt(1 - x[0] ** 2) * (1 - x[0] ** 2)),
                      2 * J * x[0] * np.sin(x[1]) / np.sqrt(1 - x[0] ** 2)]])


sol = optimize.root(func, [0.5, np.pi + 0.1], jac=jac, method='hybr')
print sol.x
print func(sol.x)