Check numerically that the constructed HO wf's are orthonormal. I.e., evaluate $\int r^2dr R_{nl}(r)R_{n'l'}(r)$. You will want to use Gaussian quadrature to discretize the integrals. Gaussian quadrature is discussed some in
. While Gauss-Laguerre quadrature is ideal for this problem, plain Gauss-Legendre quadrature, which is widely available in canned routines you can easily find via Google, is sufficient.