Normalize the solution during the integration by differentialequations.jl

I usually come to solve outwards ode like schrodinger equation under columb potential and I recently realize that the initial condition and time span affect the stability of the solution as posted as following:

using DifferentialEquations,Plots,Roots;
pyplot()

function columbode!(du,u,p,r)
    v(r)=-1/r
    du[1]=u[2]
    du[2]=(v(r)-p)*u[1]
end

f(x)=let u0=[0.0,x],rspan=(1e-10,80.0)
    prob(e)=ODEProblem(columbode!,u0,rspan,e)
    solf(e)=solve(prob(e), Tsit5(), reltol=1e-8, abstol=1e-8)
    e_eigen=find_zero(e->solf(e).u[end][1],(-1.0,-0.1))
    plot(solf(e_eigen),vars=[1])
end
display(plot(f(1e-3),f(3e-3),f(6e-3),f(1e-2),layout=(2,2)))
savefig("discourse_rescaling_ode.png")

discourse_rescaling_ode

I tested four different initial derivative and found the final solutions perform differently. Only one of them behaves normally. I wonder if there’s way to deal with stability of outwards ode. Currently, I wonder if I can normalize the intermediate sol during integration using callback or integrator interface and I have no good ideas since these operations may corrupt the interpolations stuff. Thanks for any help!

A specific question I want to ask is: is it possible to modify the entire sol.u like normalize it using callback or integrator interface?