\bigl\{ r_{ij} \bigr\})$ from Eq.(\[psi\]) is written in the form $$\begin{aligned}
 \psi_{L=0}(A; \bigl\{ r_{ij} \bigr\}) = \sum^N_{k=1} C_k r^{n_1(k)}_{23} r^{n_2(k)}_{13} r^{n_3(k)}_{12} r^{m_1(k)}_{14} 
 r^{m_2(k)}_{24} r^{m_3(k)}_{34} exp(-\alpha_{k} r_{14} -\beta_{k} r_{24} -\gamma_{k} r_{34}) \label{semexp}\end{aligned}$$ where $\alpha_k, \beta_k, \gamma_k$ ($k = 1, 2, \ldots, N$) are the varied non-linear parameters. The presence of the varied non-linear parameters in Eq.(\[semexp\]) is the main and very important difference with the traditional Hylleraas variational expansion (see, e.g., [@Lars]). In the last case in Eq.(\[semexp\]) we have $\alpha_1 = \ldots = \alpha_N, \beta_1 = \ldots = \beta_N$ and $\gamma_1 = 
\ldots = \gamma_N$. Note that all matrix elements of the Hamiltonian, Eq.(\[Hamil\]), and overlap matrix needed in computations with the use of the semi-exponential basis, Eq.(\[semexp\]), contain the same three-electron integrals which arise for the usual Hylleraas expansion (for more detail, see discussion in [@Fro2010]). In other words, numerical calculation of all matrix elements with semi-exponential functions is not a more difficult task than for the traditional Hylleraas radial functions. This also simplifies numerical computations of the bound state properties (i.e. expectation values) in the semi-exponential basis set. In particular, our algorithms used in calculations of all required matrix elements is based on the old Perkins formula for 