g[ϵ_, ξ0_, ρ_, f_]:=Evaluate[ψ[ξ0]  /. 
                  NDSolve[{-ψ''[ξ] - (ρ^2) 
                            f[ξ] ψ[ξ] == ϵ ψ[ξ], ψ[-ξ0] == 
                        1, ψ'[-ξ0] == (Sqrt[Abs])[ϵ] ψ[-ξ0]}, ψ, {ξ, -ξ0, ξ0}]] [[1]]

pl[ξ0_, ρ_, f_] := 
            Plot[Sign[
                g[ϵ, ξ0, ρ, 
                  f]], {ϵ, -ρ^2, 0}, Axes -> False, 
              Frame -> True, PlotRange -> {-1.1, 1.1}]

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

pl[15,  .3, f]

pl[15,  .5, f]

pl[15, 1, f]

pl[15, 2, f]

pl[15, 4, f]

pl[4, 10, f]

pl[5, 20, f]

g[x_] = UnitStep[1 - x^2]

pl[1.1, 20, g]