Analytical-numerical method for calculating the Thomas-Fermi potential (Q2304183)

A new method of calculating the solution of the singular nonlinear second-order boundary value problem \begin{align*} &\frac{\mathrm{d}^2}{\mathrm{d}r^2}\Psi(r)-\frac{1}{\sqrt r}\Psi^{3/2}(r)=0,\qquad r\in(0,R)\tag{1}\\ &\lim_{r\to0^+}\Psi(r)=Z,\quad\lim_{r\to R^-}\Psi(r)=0\tag{2} \end{align*} is considered. Problem (1)--(2) can be explained by modeling the spherically symmetric Thomas-Fermi problem of a heavy atom or ion at zero absolute temperature, where $r$ physically means the distance to a point nucleus with the charge $Z$, $R\in(0,\infty)$ is the radius of the ion ($R=\infty$ for an atom), and the shielding coefficient $\Psi(r)$ determines the intra-atomic Coulomb potential by the formula $\Phi(r)=\Psi(r)/r$. Note that the existence of a classical solution of Problem. (1)--(2) is known [\textit{G. Sansone}, Equazioni Differenziali nel Campo Reale. Bologna: Nicola Zanichelli (1948)] and it is uniquely determined by the pair $(R,Z)$ for $R\in(0,\infty]$. The advantage of the presented method is the high accuracy of the result achieved with relatively small computing costs. The developed method enables to obtain a solution and its derivative with a predetermined accuracy at every point of the ray $r\in[0,\infty)$ for $R=\infty$. The computational algorithm accepts at the input, instead of a pair of values $(R,Z)$, the pair $(\mu,Z)$, where $\mu\in(0,1)$ is a numerical parameter and defining the solution uniquely up to a scaling coordinate transformation for $R<\infty$. The corresponding value of $\mu$ can be established approximately via the data $R$, $Z$ and interpolation. The presented computational method uses new parametric representations of solutions of Problem (1)--(2) containing, in addition to elementary functions, also solutions to some first-order differential equations whose right-hand side is rational both in the argument and in the unknown function. These auxiliary functions admit an efficient calculation with high accuracy, since they retain analyticity along the entire interval of the parameter change, including its ends. The paper presents numerical calculations that demonstrate the capabilities of the proposed analytical-numerical methods. A detailed presentation of the results obtained here with proofs is given in [the author, Comput. Math. Math. Phys. 59, No. 8, 1292--1313 (2019; Zbl 1430.34031); translation from Zh. Vychisl. Mat. Mat. Fiz. 59, No. 8, 1358--1380 (2019); ibid. Phys. 58, No. 12, 1948--1966 (2018; Zbl 1418.34032); translation from Zh. Vychisl. Mat. Mat. Fiz. 58, No. 12, 2026--2047 (2018); ibid. 58, No. 2, 230--237 (2018; Zbl 1392.35066); translation from Zh. Vychisl. Mat. Mat. Fiz. 58, No. 2, 244--252 (2018)].