Solve a pair of coupled nonlinear equations within certain limits This answer to this question works only for situations in which the desired solution to the coupled functions is not restricted to a certain range.

But what if, for example, we wanted a solution such that 0 &lt; x &lt; 10 and 0 &lt; y &lt; 10?

There are functions within scipy.optimize that find roots to a function within a given interval (e.g., brentq), but these work only for functions of one variable.

Why does scipy fall short of providing a root solver that works for multi-variable functions within specific ranges?  How might such a solver be implemented?
As a workaround, you could minimize another function that includes both the objective and the constraints, then check if sol.fun is (numerically) equal to zero.

from scipy.optimize import minimize
import numpy as np

f = lambda x: np.sin(x).sum() #the function to find roots of
L = np.array([-1,-1]) #lower bound for each coordinate
U = np.array([1, 1])  #upper bound for each coordinate
g = lambda x: f(x) **2 +  max(0, (L - x).max(), (x - U).max())
sol = minimize(g, [0.5,0.5])


Also, scipy.optimize seems to have some optimisers that support rectangular bounds, i.e. differential_evolution (since version 0.15.0).