Normalised (Subscript[ψ, n])(ξ)for simple harmonic oscillator.

MEnO[n_,a0_,a1_,lastj_]:=Module[{},
K=2n+1;
a[1]=a1;a[0]=a0;

For[j=0,j≤lastj,j++,
a[j+2]=((2j+1)-(K))/((j+1)(j+2))a[j];
];(*end for *)

If[OddQ[n],
For[j=0,j≤lastj,j+=2,
a[j]=0;
];(*end for*)
];(*end if j*)

If[EvenQ[n],
For[j=1,j≤lastj,j+=2,
a[j]=0;
];(*end for *)
];(*end if *)
];

h[xi_,n_]:=Module[{nlast=n},
a0=1;a1=1;
jlast=nlast;
MEnO[nlast,a0,a1,jlast];
Sum[a[i]xi^i,{i,0,nlast}]
];

psi[xi_,n_]:=h[xi,n]*Exp[-(xi^2)/2.0];
(*Module[{nlast=n},
a0=1;a1=1;
jlast=nlast;
MEnO[nlast,a0,a1,jlast];
Exp[-(xi^2)/2.0]*Sum[a[i]xi^i,{i,0,nlast}]
];*)

nlast=10;

For[n=0,n≤nlast,n++,
A[n]=1.0/Sqrt[NIntegrate[psi[xi,n]^2,{xi,-Infinity,Infinity}]];
];

Npsi[xi_,n_]:=A[n]*psi[xi,n];

Manipulate[Plot[h[xi,n],{xi,-10,10},PlotRange->{-20,20},AxesLabel->{"ξ","(Subscript[h, n])(ξ)"},
PlotLabel->{"n="<>ToString[n]<>"; (Subscript[h, n])(ξ)=",h[ξ,n]}
],{n,0,nlast,1}]

Manipulate[Plot[(Npsi[xi,n])^2,{xi,-10,10},PlotRange->{-1,1},AxesLabel->{"ξ","\[LeftBracketingBar](Subscript[ψ, n])(ξ)(\[RightBracketingBar]^2)"},
PlotLabel->{"n="<>ToString[n]<>"; (Subscript[ψ, n])(ξ)=",psi[ξ,n]}
],{n,0,nlast,1}]

Print["Limit[h[xi,0],{xi->Infinity}]=",Limit[h[xi,10],{xi->Infinity}]];

For[n=0,n≤nlast,n++,
Print["n=",n];
Print["(Subscript[A, n])=",A[n]];
Print["(Subscript[h, n])(ξ)=",h[ξ,n]];
Print["(Subscript[ψ, n])(ξ)=",Npsi[ξ,n]];
Print[" "];
];

-(((d^2) ψ x)/(dx^2))+(1/2) ((mω)^2) (x^2) ψ x=Eψ x
The following code shows if n is not integer ,or equivalently, E does not take on a value E=(n+1/2)ℏω,
, the solution will diverge.

ℏ=1;m=1;ω=0.001;EE[n_]:=(n+0.5)ω*ℏ;
Solodd[n_]:=Module[{},
s=NDSolve[{-(ℏ/(2m))ψ''[x]+(1/2)m (ω^2)*(x^2)*ψ[x] ==EE[n]*ψ[x],ψ[0]==0,ψ'[0]==-1},ψ,{x,-1000,1000}];
dum=Evaluate[ψ[x]/.s];
Plot[{dum/500},{x,-1.0xmax,1.0xmax},PlotStyle->Automatic,PlotRange->{-0.2,0.2},PlotLabel->"n="<>ToString[1.0n]<>"; E=(n+1/2)ℏω"
,AxesLabel->{"x","(Subscript[ψ, n])(x)"}]
];

Soleven[n_]:=Module[{},
s=NDSolve[{-(ℏ/(2m))ψ''[x]+(1/2)m (ω^2)*(x^2)*ψ[x] ==EE[n]*ψ[x],ψ[0]==1,ψ'[0]==0},ψ,{x,-1000,1000}];
dum=Evaluate[ψ[x]/.s];
Plot[{dum/10},{x,-1.0xmax,1.0xmax},PlotStyle->Automatic,PlotRange->{-0.2,0.2},AxesLabel->{"x","(Subscript[ψ, n])(x)"},
PlotLabel->"n="<>ToString[1.0n]<>"; E=(n+1/2)ℏω"
]
];

Sol[n_]:=If[EvenQ[n],Soleven[n],Solodd[n]];
xmax=180;
Print[Manipulate[Sol[frame],{frame,0,10,1/2}]]