Rayleigh-Ritz method in J. Thijssen, Computational Physics 2nd edition, (Cambridge University Press, 2007)

Rayleigh-Ritz.jl

using LinearAlgebra
α = [13.00773, 1.962079, 0.444529, 0.1219492]
S = [(π/(α[i]+α[j]))^(3/2) for i in keys(α), j in keys(α)]
H = [3*π^(3/2)*α[i]*α[j]/(α[i]+α[j])^(5/2) - 2*π/(α[i]+α[j]) for i in keys(α), j in keys(α)] 
E, C = eigen(H, S)

julia> using LinearAlgebra

julia> α = [13.00773, 1.962079, 0.444529, 0.1219492]
4-element Vector{Float64}:
 13.00773
  1.962079
  0.444529
  0.1219492

julia> S = [(π/(α[i]+α[j]))^(3/2) for i in keys(α), j in keys(α)]
4×4 Matrix{Float64}:
 0.0419641  0.0961392   0.112858   0.117043
 0.0961392  0.716317    1.49148    1.85084
 0.112858   1.49148     6.64247   13.0602
 0.117043   1.85084    13.0602    46.2287

julia> H = [3*π^(3/2)*α[i]*α[j]/(α[i]+α[j])^(5/2) - 2*π/(α[i]+α[j]) for i in keys(α), j in keys(α)]
4×4 Matrix{Float64}:
  0.577268    0.0720025  -0.32154    -0.436126
  0.0720025   0.50705    -0.989185   -2.37742
 -0.32154    -0.989185   -2.63808    -7.34222
 -0.436126   -2.37742    -7.34222   -17.3052

julia> E, C = eigen(H, S)
GeneralizedEigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
4-element Vector{Float64}:
 -0.4992784056674856
  0.11321392045798631
  2.592299571959811
 21.144365190122553
vectors:
4×4 Matrix{Float64}:
 -0.517825   0.240729  -0.00593851   1.0
 -0.878387   0.163899   1.0         -0.201492
 -1.0        1.0       -0.360593     0.0367844
 -0.397122  -0.414971   0.0560347   -0.00500071

J. Thijssen, Computational Physics 2nd edition, (Cambridge University Press, 2007)

$$ \begin{aligned} \phi_1(r) &= \exp(-\alpha_1 ~r^2) = \exp(-13.00773 ~r^2), \\ \phi_2(r) &= \exp(-\alpha_2 ~r^2) = \exp(-1.962079 ~r^2), \\ \phi_3(r) &= \exp(-\alpha_3 ~r^2) = \exp(-0.444529 ~r^2), \\ \phi_4(r) &= \exp(-\alpha_4 ~r^2) = \exp(-0.1219492 ~r^2). \end{aligned} $$

$$ \begin{aligned} S_{ij} = \langle \phi_{i} | \phi_{j} \rangle &= \int \phi_{i}^*(r) \phi_{j}(r) \mathrm{d} \pmb{r} \\ &= \underline{\left( \frac{\pi}{\alpha_i + \alpha_j} \right)^{3/2}} \end{aligned} $$

$$ \begin{aligned} V_{ij} = \langle \phi_{i} | \hat{V} | \phi_{j} \rangle &= \iiint \phi_{i}^*(r) \left[ -\frac{1}{r} \right] \phi_{j}(r) ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \underline{- \frac{2\pi}{\alpha_i + \alpha_j}} \end{aligned} $$

$$ \begin{aligned} T_{ij} = \langle \phi_{i} | \hat{T} | \phi_{j} \rangle &= \iiint \mathrm{e}^{-\alpha_i r^2} \left[ -\frac{1}{2} \nabla^2 \right] \mathrm{e}^{-\alpha_j r^2} ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \underline{ \frac{3 \alpha_i \alpha_j \pi^{\frac{3}{2}}}{(\alpha_i + \alpha_j)^{\frac{5}{2}}} } \end{aligned} $$

$$ \begin{aligned} H_{ij} &= \langle \phi_{i} | \hat{H} | \phi_{j} \rangle \\ &= \langle \phi_{i} | \hat{T} + \hat{V} | \phi_{j} \rangle \\ &= \langle \phi_{i} | \hat{T} | \phi_{j} \rangle + \langle \phi_{i} | \hat{V} | \phi_{j} \rangle \\ &= \underline{T_{ij} + V_{ij}} \end{aligned} $$