## 3 Solution Procedure
### Scheme for Numerical Solution for the System of Equations for MO
The solution of Equation (21) must satisfy the smoothness requirement corresponding to the equality of the first-order derivatives of a function obtained by sweeping the values \"on the left\" and \"on the right\" at the crosslinking point \(Z\):  
\[\left(P_{Z}^{L}\right)^{\prime}=\left(P_{Z}^{R}\right)^{\prime} \tag{27}\]  
This condition leads to a system of linear inhomogeneous equations with respect to the value of the partial harmonic at the crosslinking point \(P_{Z}\):  
\[\big{(}\Big{[}K_{Z}^{L}\Big{]}^{-1}-K_{Z+1}^{R}\big{)}P_{Z}=M_{Z+1}^{R}+\Big{[} K_{Z}^{L}\Big{]}^{-1}M_{Z}^{L}. \tag{28}\]  
The algorithm for implementing the numerical procedure of solving the system of Equation (21) consists in determining the eigenvalues \(\varepsilon_{n\gamma}\) and corresponding MO of the \(\gamma\)-type symmetry. The solution of the homogeneous system of equations obtained from Equation (21) by zeroing out the nonlocal summand \(X=0\) is performed at the stage of initialization of the iterative procedure using the initial approximation (see Section 3.2) for MO. This is, in fact, the solution of the Hartree equation, which does not take into account the exchange terms. When the numerical procedure is implemented, the homogeneous Equation (28) is solved, which right-hand side is equal to zero. Moreover, the existence of non-zero solutions is determined by the criterion:  
\[\mathrm{det}\|\Big{[}K_{Z}^{L}\Big{]}^{-1}-K_{Z+1}^{R}\|=0. \tag{29}\]  
The selection of solutions on the energy scale until reaching the corresponding MO \(n\gamma\). For example, the \(\sigma\) symmetry states for diatomic molecules can be selected by the principal quantum number. In this case, the search for highly excited states of the discrete spectrum involves the sequential determination of all lower one-electron energies. The implemented procedure for finding solutions accounts the obtained spectrum of state energies and uses them as \"reference points\" on the one-electron energy scale \(\varepsilon_{n\gamma}\) of MO within the range of \([-E,0]\) at subsequent iterations of the computational procedure. This approach ensures time saving for calculations and increases the efficiency of the algorithm for finding solutions of Equation (21). The final set of radial parts of the partial harmonics \(P_{I}^{n\gamma}(r)\ c\ l\leq l_{0}\) is found for MO when solving Equation (21). Following the procedure given in Section 2.2, the partial harmonics with orbital angular momenta \(l>l_{0}\) are determined by the superposition of ligand functions. The participation coefficients \(\Lambda_{\mu\alpha}^{n\gamma L}\) of the ligand wave functions \(\Psi_{ul_{0}}^{n\gamma L}(r)\) are determined by the solution of the system of algebraic Equation (14). Here, the functions \(\Psi_{\mu l_{0}}^{n\gamma L}(r)\) are preformed according to the numerical procedure given in Section 2.3 on the basis of the wave functions of electron shells of atoms, obtained by solving the Hartree-Fock equations. Equation (21) has a unique solution for any value of one-electron energy \(\varepsilon_{n\gamma}\) when accounting nonlocal terms \(X\). Therefore, the selection of solutions is based on the overlap integral, which is obtained according to the expression:  
\[\langle\Phi_{n\gamma}|\Phi_{n\gamma}\rangle=\sum_{l}\int P_{n\gamma l}^{2}(r)dr \tag{30}\]  
The criterion for selecting solutions for MO at the intermediate iteration step of the computational procedure is the correspondence of the overlap integral of the wave function (Equation (30)) to the value obtained at the first iteration of the matching process.