The energy difference between $H$ atom and $H^+$ ion in Thomas Fermi Theory

Question

The famous Thomas Fermi theory says that we can express the total energy of a atomic system with $N$ electrons orbiting around a charged nucleus of atomic number $Z$ in terms of radial density of electrons $\rho(r)$ as,
$E_{Total}=K\int_{\mathbb{R}^3} \rho(r)^{5/3} dr + \frac{e^2}{2}\int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \frac{\rho(r) \rho (r')}{|r-r'|} dr dr'-Ze^2\int_{\mathbb{R}^3} \frac{\rho(r)}{|r|} dr $
Where the 1st term represents the kinetic energy of the electrons ($K$ is a constant), the 2nd term represents the electron-electron repulsion energy and the last term is the energy due to the attraction due to the positively charged nucleus and the electrons.
Now my question is how can I find out the energy difference between the neutral Hydrogen atom ($H$) and the positive ion of Hydrogen ($H^+$) and similar for $H$ and $H^-$ ion?
I have tried out to minimize the energy using Lagrange multiplier and arrived at $\mu=\frac{5}{3}K \rho(r)^{5/3}+e^2\int_{\mathbb{R}^3}\frac{\rho(r')}{|r-r'|}dr'-Z\frac{e^2}{|r|}$. My confusion here is how to find out the density of electrons $\rho(r)$ in case of Hydrogen and any atom to derive the total energy of the system.

Edit:
Defining $V(r)=-\frac{Ze}{|r|}$ and $U(r)=e\int_{\mathbb{R}^3} dr'\frac{\rho(r')}{|r-r'|}$, we can write the total potential as $V_{tot}=V(r)+U(r)$. Now, $$\nabla^2 V_{tot}=-4\pi \rho(r)$$ Using sperical symmetry we write $\nabla^2=\frac{1}{r^2}\frac{\partial}{\partial r}\big(r^2 \frac{\partial}{\partial r}\big)$. And we have already arrive at $\rho(r)=\big(\frac{3}{5K}\big)^{3/2}\big[\mu-eV_{tot}(r) \big]^{3/2}$. Finally using $\mu=0$ and $\rho(r)=-\frac{1}{4\pi}\nabla^2V_{tot}(r)$ we get the following differential equation, $$\frac{d^2 V_{tot}(r)}{dr^2}+\frac{2}{r}\frac{dV_{tot}(r)}{dr}+QV_{tot}(r)^{3/2}=0$$ here $Q$ is a constant.
If we assume $V(r)\sim r^{-\alpha}$, we get $\alpha =4$ and $\rho(r) \sim r^{-6}$
Now the questions are:
1. Is this expression of $\rho(r)$ correct?
2. Do we have any better method to solve the differential equation?

Answer 1

In a positive hydrogen ion, there are no electrons, so the electron density is zero.

In the other cases, you seem to be assuming that your problem is spherically symmetric, which implies that everything is in the ground state. In that case, you can safely approximate the electron charge density with $e|\psi_{100}|^2$, where the latter term is the ground state of the hydrogen atom,

$\psi_{100}=\frac{1}{\sqrt{\pi a_0^3}}e^{-r/a_0}$

where $a_0$ is the Bohr radius.

For the negative ion, there are two electrons now, so the situation is more complicated. There are actually no analytic solutions for the wave function, but by variation of parameters one can approximate the wave function as (before normalization)

$\psi(r_1,r_2)=e^{-(ar_1+br_2)}+e^{-(br_1+ar_2)}$

where the best parameter values (i.e. those that minimize the energy of the system) are $a=1.04,b=0.28$. These results came from https://arxiv.org/pdf/0907.2614.pdf, and at this point the electron density can be extracted from the wavefunction in a similar manner, averaging over the position of the two electrons.