Solution (source code)

= Solution

Enlarge the trial space while retaining the optimized Gaussian. For example, use
$$
\psi(x)=e^{-\alpha x^2/2}(c_0+c_1x^2+c_2x^4+\cdots+c_Nx^{2N}),\qquad N\geq2,
$$
normalize it and minimize the <Rayleigh quotient> over its coefficients, optionally optimizing $\alpha$ too. This even trial space respects the even ground-state parity of the confining symmetric potential and adds independent shape freedom beyond a Gaussian width.

The original Gaussian remains included by setting the later coefficients to zero, so the optimized upper bound cannot increase. Finite-basis optimization is the <Rayleigh-Ritz method>, with matrix entries given by the Hamiltonian and overlap integrals. Increasing the basis size systematically improves the available approximation.