Numerical solution of the Schroedinger equation for the harmonic oscillator.

a=-10;
b=10;
c=-5;
d=5;

ϵ is the energy in units of (hbar ω/2)

ϵ=0.999999964931055;

Even solutions

Plot[
Evaluate[f[u]/.NDSolve[{f''[u] - (u^2 - ϵ)f[u]==0,f[0]==1,f'[0]==0},f[u],{u,-10,10},MaxSteps->10000
]],{u,a,b},PlotRange->{c,d}
]

ϵ is the energy in units of (hbar ω/2)

ϵ=2.999999991504038;

Odd solutions

Plot[
Evaluate[f[u]/.NDSolve[{f''[u] - (u^2 - ϵ)f[u]==0,f[0]==0,f'[0]==1},f[u],{u,-10,10},MaxSteps->10000
]],{u,a,b},PlotRange->{c,d}
]