>>> 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])