Solve a pair of coupled nonlinear equations within certain limits
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).