Discretized Harmonic Oscillator H = -(p^2)/(2 m) + 1/2 m ω^2 (x^2)⟶ -1/2 ψ''+(1/2)(s^2)ψ== ϵψ

discretize on grid of dx=0.05, |s|<0.5

ss=Range[-5,5,.05];

in this basis, potential energy matrix is diagonal

V=DiagonalMatrix[(1/2)(ss^2)];

V[[10,;;10]]//MatrixForm

approximating ψ''(s)=(ψ(s-ds)+ψ(s+ds)-2ψ(s))/((ds^2))

KE=((-1)/(2 (((.05))^2)))(-2 DiagonalMatrix[1+0 Range[201]]+ DiagonalMatrix[1+0 Range[200],1]+ DiagonalMatrix[1+0 Range[200],-1]);

H[[;;10,;;10]]//MatrixForm

H=KE+V;

{eval,kets}=Eigensystem[H];

eval[[-10;;-1]]

note excellent agreement of eigenvalues with analytical theory

plot of first 3 states

Show[Table[ListPlot[{ss,kets[[m]]}\[Transpose],Joined->True,PlotRange->All,PlotStyle->Hue[-m/3]],{m,-3,-1}]]