System of nonlinear equations:
$$E_i=\epsilon_i+\sum_{j\neq i}^N \left(\frac{1}{1+\exp(E_j/T)}-\frac{1}{2}\right)\frac{e^2}{r_ij} \tag 1$$
where $T=0.05$, $r_{i,j}$ is given symmetric matrix with trace $0$ and $\epsilon_i$ is random vector with uniform distribution in the range $[-0.5, 0.5]$
How do I find the metastable solutions for that system of equations efficiently i.e. without using newton-raphson or some similar methods (let's say $i$ and $j$ can go from $1$ to $10000$)?
Proposed solution: (don't have to read if you don't have power for it)
Setting eq. $(1)$ in a units of $\frac{e^2}{r_{nn}}$ where $r_{nn}$ is some constant we get
$$E_i=\epsilon_i+\sum_{j\neq i}^N n_j\frac{r_{nn}}{r_ij}\tag2$$
$$n_j=\frac{1}{1+\exp(E_j/T)}-\frac{1}{2}\tag3$$
Starting with random uniform distribution of $n_j$ in the range $[-0.5, 0.5]$ and finding $E_i$ though eq. $(2)$ and then input the resulted $E_i$ in $(3)$ to find new $n_j$ and feed it back in $(2)$ again to find new $E_i$. After alot of iterations I think I suppose to get some metastable solution under some specific configuration of $E_i$ and $n_i$ which ofcourse will stongly be dependet at the initial values of $\epsilon_i$ and $n_j$.
The reason I'm doing this and not using newton-raphson (or similar methods like fsolve function in matlab) is because $E_i$ and $n_j$ are huge vectors and the computation time will be extremely long.
*The problem that I have with this method is that I cant find a way to preserve a very importent property, $\sum_jn_j=0$, I'm losing it after the first iteration.