In quantum mechanics we can maybe express the s-states (spherically symmetric wave functions) of the Helium atom as two wave function depending on the spherical radius variable $r$ ($0 < r < + \infty$) which is the distance from the nucleus. In this formulation of Helium type atoms we have the coupled system of ODE:
$$ \left[ - \frac{\partial_r^2}{2} - \frac{Z}{r} + \frac{1}{2} \left( \int_0^r R_2(u)^2 d u + r \int_r^{+ \infty} \frac{R_2(u)^2}{u} d u \right) \right] R_1(r) = \epsilon_1 R_1(r) $$ $$ \left[ - \frac{\partial_r^2}{2} - \frac{Z}{r} + \frac{1}{2} \left( \int_0^r R_1(u)^2 d u + r \int_r^{+ \infty} \frac{R_1(u)^2}{u} d u \right) \right] R_2(r) = \epsilon_2 R_2(r) $$
With $R_1(r)$, $R_2(r)$ the wave functions of the electron 1 and electron 2. Then $\epsilon_1$ and $\epsilon_2$ are the corresponding eigenvalue of each (their energy in physics). $R_1$ and $R_2$ can be equal for certain case (two $1s$ state as ground state, two $2s$ states, two $3s$ states and so on for spherically symmetric excited states).
To this we add boundary condition, we have only these ones:
$$ R_i(0) = 0$$ $$ R_i(+ \infty) = 0 $$
and the normalization condition defined as:
$$\int_0^{+ \infty} R_i(r)^2 d r = 1$$
with ${i = 1,2}$. I know also that we can "cheat" on boundary condition by putting $\partial_r R_i(0) = 1$ or whatever positive constant. It works for Hydrogen atom because the equation is linear so I do not know if it works here because it is non-linear.
QUESTION:
How to solve for eigenvalues $\epsilon_1$, $\epsilon_2$ and eigenfunctions $R_1(r)$, $R_2(r)$ with these boundary conditions using the shooting method ?
