Solving Schrödinger's Equation for the electronic energies of the Molecular Ion Hydrogen H2+ in the Elliptic coordinate system

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

As you have discovered, the H$_2^+$ ion is separable but it is not exactly solvable. Spheroidal coordinates allow you to separate the three-dimensional time-independent Schrödinger equation into three separate one-dimensional Schrödinger problems, one of which is trivial, but that's as far as you can go.
Numerically that's a complicated problem, because the equations are not quite as separated as you'd like: the eigenvalue $\varepsilon$ and the separation constant $A$ appear in both the $\xi$ and the $\eta$ equations, so they make a 'bispectral' problem and must be solved in tandem.