Here is a symbolic solution of your eigenvalue problem.
Define the differential equation (setting $\hbar = \omega = m_0 = 1$).
diffeq = -(1/2) \[Psi]''[x]ψ''[x] + 1/2 x^2 \[Psi][x]ψ[x] == e \[Psi][x]ψ[x]
Symbolically solve the differential equation.
soln = DSolve[diffeq, \[Psi]ψ, x][[1, 1]]
(* \[Psi]ψ -> Function[{x},
C[2] ParabolicCylinderD[1/2 (-1 - 2 e), I Sqrt[2] x] +
C[1] ParabolicCylinderD[1/2 (-1 + 2 e), Sqrt[2] x]] *)
Inspect the $x \longrightarrow +\infty$ limit by series expanding there:
\[Psi][x]ψ[x] /. soln // Series[#, {x, \[Infinity]∞, 1}] & // Simplify //
Normal // Collect[#, E^(x^2/2), Simplify] &
(* 2^(-(1/4) + e/2) E^(-(x^2/2)) Sqrt[1/x] x^
e C[1] + (1 - I) 2^(-(3/4) - e/2) E^(x^2/2) Sqrt[1/x]
x^-e C[2] (Cos[(e \[Pi]π)/2] - I Sin[(e \[Pi]π)/2]) *)
We need C[2]->0 in order to suppress the divergent exponential term.
Now inspect the $x \longrightarrow -\infty$ limit by series expanding there:
\[Psi][x]ψ[x] /. soln /. {C[2] -> 0} //
Series[#, {x, -\[Infinity]∞, 1}] & // Simplify // Normal // Expand
(* -2^(-(1/4) + e/2) E^(2 I e \[Pi]π - x^2/2) (1/x)^(1/2 - e) C[1]
- (I 2^(1/4 - e/2) E^(-I e \[Pi]π + x^2/2) Sqrt[\[Pi]]Sqrt[π] (1/x)^(1/2 + e)
C[1])/Gamma[1/2 - e] *)
We need e -> n + 1/2 (where n is a non-negative integer) to make the 1/Gamma[1/2 - e] suppress the divergent exponential term. This makes use of the fact that $\Gamma(x)$ is singular for $x = 0,-1,-2,-3,\cdots$.
Gathering it all together gives the solution
\[Psi][x]ψ[x] /. soln /. {C[2] -> 0, e -> n + 1/2}
(* C[1] ParabolicCylinderD[1/2 (-1 + 2 (1/2 + n)), Sqrt[2] x] *)
which can be normalised by choosing C[1] appropriately.
