Numerical solution of harmonic oscillator V=(1/2)(bx^2)

solve Schrodinger eqn,

SE[ϵ_]:=((-1)/2)psi''[s]+(1/2)(s^2)psi[s]==ϵ psi[s]

boundary condition deep into forbidden region

BC[s0_]:={psi[s0]==1,psi'[s0]==-s0 psi[s0]}

Numerically integrate SE from y0 in to 0, calculate values at y=0
Try different values of ϵ to either get psi’[0]=0 or psi[0]=0, depending on whether you are after symmetric or antisymmetric solution
Note:  If the potential was not symmetric, then you would also need to integrate from -y0 to 0 and the condition for a proper eigenfunction would be that the wavefunction and its derivative would be continuous at y=0.

With[{ϵ=34.5,s0=50},
(*Numerical Integration of SE from y0 to 0*)
sol=NDSolve[{SE[ϵ],BC[s0]},psi,{s,s0,0}][[1]];
(*print wavefunction and derivative at 0*)
Print[{((psi'[0])/(psi[0])),((psi[0])/(psi'[0]))}/.sol];
(*plot result*)
Plot[((psi[y])/(psi[0]))/.sol,{y,0,15},PlotRange->All]
]