>>> from scipy.optimize import fsolve
>>> def func(x):
...     return [x[0] * np.cos(x[1]) - 4,
...             x[1] * x[0] - x[1] - 5]
>>> root = fsolve(func, [1, 1])
>>> root
array([6.50409711, 0.90841421])
>>> np.isclose(func(root), [0.0, 0.0])  # func(root) should be almost 0.0.
array([ True,  True])