Electronic Energies of Molecular Ion Hydrogen $H_2^{+}$
$r_1$ is the distance between the proton $1$ and the electron.
$r_2$ is the distance between the proton $2$ and the electron.
$R$ is the distance between the two protons, fixed parameter in Born Oppenheimer approximation.
Schrödinger's Equation:
$$ \hat{H} \Psi = E_{el} \Psi $$
$$ -\frac{{\hbar}^2}{2 m_e} \Delta \Psi -\frac{e^2}{4 \pi \varepsilon_0} \left(\frac{1}{r_1}+\frac{1}{r_2}\right)\Psi = E_{el} \Psi $$
Atomic Units:
$$ a_0=\frac{4 \pi \varepsilon_0 \hbar^2}{m_e e^2} = 0.5291 *10^{-10} m $$
$$ E_h=\frac{\hbar^2}{m_e a_0^2}=27.21 eV $$
$$ \Delta \Psi + \frac{2}{a_0}\left(\frac{1}{r_1}+\frac{1}{r_2}\right)\Psi=-2 \frac{E_{el}}{E_h a_0^2} \Psi $$
We introduce the adimentional parameters $ r= \frac{R}{a_0} $ and $ \varepsilon = \frac{E_{el}}{E_h} $
Elliptic coordinate system:
$ \xi = \frac{r_1 + r_2}{R} $ , $ \eta=\frac{r_1 - r_2}{R} $ and $ \phi $ is the angle between the position vector of the electron and the $xy$ plane.
The Laplacian operator will become:
$$ \Delta =\frac{4}{R^2 (\xi^2 - \eta^2) } \left[\frac{\partial }{\partial \xi} (\xi^2 -1)\frac{\partial }{\partial \xi}+\frac{\partial }{\partial \eta} (1-\eta^2)\frac{\partial }{\partial \eta}+\frac{\xi^2 - \eta^2}{(\xi^2-1)(1-\eta^2)}\frac{\partial^2 }{\partial \phi^2}\right] $$
then we obtain:
$ \left[\frac{\partial }{\partial \xi} (\xi^2 -1)\frac{\partial }{\partial \xi}+\frac{\partial }{\partial \eta} (1-\eta^2)\frac{\partial }{\partial \eta}+\frac{\xi^2 - \eta^2}{(\xi^2-1)(1-\eta^2)}\frac{\partial^2 }{\partial \phi^2}\right] \Psi + 2r\xi \Psi = - \frac{1}{2} r^2\varepsilon (\xi^2 - \eta^2)\Psi $
We can separate the variables using $ \Psi = \Xi_{\xi} H_{\eta} \Phi_{\phi} $ so we can obtain three indipendent differential equations:
$$ \frac{\partial }{\partial \xi} (\xi^2 -1)\frac{\partial \Xi }{\partial \xi} + \left(2r\xi+A+\frac{1}{2}r^2 \varepsilon \xi^2 - \frac{\Lambda^2}{\xi^2-1}\right)\Xi=0 $$
$$ \frac{\partial }{\partial \eta} (1-\eta^2)\frac{\partial H }{\partial \eta} + \left(-A-\frac{1}{2}r^2 \varepsilon \eta^2 - \frac{\Lambda^2}{1-\eta^2}\right)H=0 $$
$$ \frac{\partial^2 \Phi }{\partial \phi^2} = - \Lambda^2 \Phi $$
where A is a separation parameter and $ \Lambda \in \mathbb{N} $
focusing on the first one we can divide both bembers for $ \Xi $ and do this substitution: $ f_{\xi} =- \frac{1}{\Xi} \frac{\partial \Xi}{\partial \xi} $ obtaining this new differential equation:
$$ 2 \xi f + (\xi^2-1)(f'+f^2)-\left(2r\xi+A+\frac{1}{2}r^2 \varepsilon \xi^2 - \frac{\Lambda^2}{\xi^2-1}\right)=0 $$
and from now on I don't know how to move. Is there someone who knows how to solve numerically this equation? Thank you very much