With scalar functions it is easy to construct a mixed Newton-Raphson-bisection algorithm so that the solution always stays inside the given bounds in which it is bracketed. What if I have a vector function which I want to finds its roots (each root depending on the other roots), but I know that the true roots have the same sign as my initial guess, and that it must stay that way in every step of the iteration procedure in order that the solution doesn't become irretrievable. Can I use this information to ensure the convergence of my solution?
To be more clear, I have a vector $\pmb{f}$ of functions, which are all the same except evaluated at different positions. Each of these component functions depend on a root vector $\pmb{x}$ comprised of all the roots at these positions. Symbolically, $\pmb{f}(\pmb{x})$ looks like $$f_i(\pmb{x})=\int_{y_i-a}^{y_i}g(y^{\prime})x(y^{\prime})dy^{\prime}+h(y_i)$$ where $x=x(y)$ is a function of position. Equivalently (?), $\pmb{x}$ can be considered as a vector solution where each of the different discrete values $y$ takes gives a component $x_i=x(y_i)$.
Now I know this is an integral equation system, but mine is actually much more complex and I'd like to solve it with Newton-Raphson if possible. The problem is that with NR alone, the solution tends to wander off in a region where some of the roots change sign, and this makes sure the solution gets forever lost.
Just to make sure, the way I proceed using only Newton-Raphson to find the roots of a vector function $\pmb{f}(\pmb{x})$ is by computing $$dx_i(\pmb{x})=\frac{f_i(\pmb{x})}{f_i^{\prime}(\pmb{x})}$$ where $$f_i^{\prime}(\pmb{x})=\frac{\partial f_i(x)}{\partial x}\Bigr|_{\pmb{x}}$$ where the boldface was dropped to show the differentiation is with respect to a variable, and then applying $$\pmb{x}_{n+1}=\pmb{x}_n-\pmb{dx}_n$$ until convergence.
I tried to use some kind of bisection method (which I know doesn't work for multidimensional problems) but I can't seem to be able to bracket a root. I asked that ${x_1}_i={x_i}_0/10$ and ${x_2}_i={x_i}_0\times 10$, calculated $\pmb{f}_1=\pmb{f}(\pmb{x}_1)$, for example. But it doesn't seem to me that such method can be generalized for a vector function.
In a way, I guess I could try a multivariate Newton-Raphson method, but usually in these cases, the vector $\pmb{x}$ is a vector of distinct variables. In my case, it's a single variable but which depends on position, hence is numerically an array. Doing partial derivatives would mean doing a derivative of that variable only at a specific position, and since I have to integrate over all positions (see below for the form of the $f_i(\pmb{x})$), it would be pretty awkward I think. Or maybe I am mistaken and it has nothing to do with Newton-Raphson. At least I know the method I explained worked for a relatively simple case, as it gave the same result as another root-finding scheme.
Is there a way to solve such a problem? I guess that if I try to solve the integral equation problem I am not guaranteed that the solution won't change sign in the process either?
