## 3 Solution Procedure  
### Scheme for Numerical Solution for the System of Equations for MO  
The numerical procedure of MO finding is reduced to solving a system of coupled differential Equation (2) for the radial parts of partial harmonics \(P_{l}^{n_{T}}(r)\) of Equation (1). The Equation (2) in general can be represented in a matrix form:  
\[P^{n}=\Omega P+X, \tag{21}\]  
where \(\Omega,X\) are two-dimensional and one-dimensional matrices respectively, each element of which is a real function specified in the coordinate grid; \(P\) is a one-dimensional matrix of the MO partial harmonics. Three-point difference Numerov scheme is applied to solve the equation [30; 31]:  
\[\mathrm{A}_{n+1}P_{n+1}-\mathrm{B}_{n}P_{n}+\mathrm{A}_{n-1}P_{n-1}=J_{n}+O \left(h^{6}\right), \tag{22}\]  
where \(\mathrm{A}_{n}=\mathrm{I}-\frac{h^{2}}{\mathrm{I}2}\Omega_{n};\mathrm{B}_{n}= 2\mathrm{I}-\frac{10\ \mathrm{h}^{2}}{\mathrm{I}2}\Omega_{n},J_{n}=\frac{h^{2}}{ \mathrm{I}2}(X_{n+1}+10X_{n}+X_{n-1});\ \mathrm{I}\) is the identity matrix; and \(n\) is the number of point at the coordinate grid. The exponential decay of the partial harmonics of MO as \(\mathbf{r}\rightarrow\mathbf{0}\) (\(\mathbf{r}\rightarrow\infty\)) allows to reduce the Cauchy problem for a second order differential equation to the boundary-value problem, using the [PERSON] algorithm [32] for the solution. This algorithm connects the values of a function at neighboring points in the following way:  
\[P_{n}=K_{n}P_{n+1}+\mathrm{M}_{n}. \tag{23}\]  
Substituting Equation (17) into Equation (16) enables obtaining the following recurrence relations:  
\[K_{n}=D_{n}^{-1}\mathrm{A}_{n+1};\,M_{n}=D_{n}^{-1}(\mathrm{A}_{n+1}\mathrm{M }_{n-1}-G_{n});\,D_{n}=\mathrm{B}_{n}-\mathrm{A}_{n-1}K_{n-1}. \tag{24}\]  
Finding a solution is carried out by independent calculation of coefficients \(K_{n}\) and \(M_{n}\) sweeping the values in the range of \(n=1\dots Z\), <<on the left>> from the selected crosslinking point \(Z\) and <<on the right>> by reversing the values \(n=N\dots Z\). Thus, an independent set of coefficient valuescan be obtained <<on the left>>\(K_{n}^{L},\ M_{n}^{L}\) and <<on the right>>\(K_{n}^{R},\ M_{n}^{R}\) from the crosslinking point \(Z\) of the solutions. Boundary conditions are used in order to obtain take-off ratios of recurrence relationships (Equation (24)). The conditions are imposed on the radial parts of partial harmonics of MO and are determined by their exponential attenuation if \(r\to 0\ (r\rightarrow\infty)\). Application of such conditions to Equation (24) results in matrix equations  
To obtain the starting coefficients of the recurrence relations (Equation (24)), the boundary conditions are applied on the radial parts of the MO partial harmonics, which are due to their exponential attenuation as r \(\to 0\ (r\rightarrow\infty)\). Applying these conditions to (24) leads to matrix equations in the form:  
\[K_{2}^{L}=\Big{(}\mathrm{B}_{2}-\mathrm{A}_{1}K_{2}^{L}\Big{)}\mathrm{A}_{3} \tag{25}\]  
\[K_{N-1}^{R}=\Big{(}\mathrm{B}_{N-1}-\mathrm{A}_{N}K_{N-1}^{R}\Big{)}\mathrm{A} _{N-2}. \tag{26}\]  
The iterative solution of Equations (25) and (26) allow restoring the coefficient values \(K_{n}^{L},\ M_{n}^{L}\), and \(K_{n}^{R},\ M_{n}^{R}\) at all points of the coordinate grid.