In the course of solving the time independent quantum harmonic oscillation Schrodinger equation $$ \Psi^{ \prime \prime} (y) +(2 \epsilon -2y^2) \Psi (y)$$
When we try ansartz $\Psi = u(y) e^{-\frac{y^2}{2}} $, we get Hermite Differential equation $$ u^{\prime \prime} - 2y u^{\prime} + (2 \epsilon -1) u $$
which admit two linearly independent fundamental power series solution, one consists of only even powers and the other consists of only odd powers.
My question is how to derive the quantization condition $$2 \epsilon - 1 = 2n $$
rigorously from the condition for 'Physical Hilbert Space' $$ \int_{-\infty}^{\infty} | u(y) |^2 e^{-y^2} dy < \infty $$
In the text books and some lecture notes I have access, they all recourse to some kind of approximation analysis that the Hermite function behaves like $e^{y^2 } $ in the limit $y \to \infty$ unless they are truncated by the quantization of $\epsilon$ but for me none of them seem mathematically exhausative. They just seems.... heuristic.
Thanks in advance.