When I use the newton's method or hybrid solver in the GSL package to deal with 1-D or multidimensional root solving problems, the code frequently crashes when the solver requests function value outside its domain of definition. For example, to solve $f(x)=0$ where $f(x)$ is only defined at $x\geq0$, even when I start with an initial guess very close to the root, the solver may still ask to evaluate the $f(x)$ at a negative $x$. Bisection method may be helpful in the 1-D problem, but it won't work for a multidimensional problem.
I typically try to solve the problem by arbitrarily define the $f(x)$ in the whole domain. But in some situation, especially in the case of complicated multidimensional root solving, I feel it's hard to extend the domain of the function and make sure it roughly maintain the general trend of the original function.
So I'm wondering if there is a method to restrict the region that solver would evaluate the function.