This note book is for numerical solutions of even eigenstates of symmetric potentials .   It is set up for a Gaussian but htis can be altered by changing the form of f[x_].

f[x_] = (E^(-x^2))

ρ = 10
ξ0 = 2
ϵ = -91
w = 
        NDSolve[{-ψ''[ξ] - (ρ^2) 
                  f[ξ] ψ[ξ] == ϵ ψ[ξ], ψ[0] == 1, ψ'[0] == 0}, ψ, {ξ, 
            0, ξ0}]
Plot[Evaluate[ψ[ξ] /.  w], {ξ, 0, ξ0}, 
        PlotRange -> {0, 5}]