Weird Schrodinger solution by DifferentialEquations.jl versus time interval

Suppose we want to solve a simple schrodinger equation under columb potential,namely,
\Psi''=(V-E)\Psi and V(r)=-\frac{1}{r},
The wavefunction u(r) of ground state is: \sim a\times r\cdot e^{-br} which should be the exact solution no matter how far the interval.But the numerical solution turns to be weired. For very large r_{max}, the wavefunction goes to large again for a certain point.

The codes goes as following:

 using DifferentialEquations,Plots,Roots
    function columbode!(du,u,p,r)
        v(r)=-1/r
        du[1]=u[2]
        du[2]=(v(r)-p)*u[1]
    end
    
    function plotf(rmax)
        u0=[0.0,1.0]
        rspan=(1e-10,rmax)
        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
    
    plot(plotf(50.0),plotf(60.0),plotf(70.0),plotf(80.0),layout=(2,2))

The result is as following:
捕获

Any idea on how to solve this problem?By the way, I don’t want to use callback to set the wavefunction to artificial zero after I need to find the root of wavefunction with respect to trial energy. Thanks!

The outwards equation is very unstable (if you mess up even a little bit on E or u0, the solution will blow up). Either tighten the tolerance further, use a small interval or possibly use a more stable variant of this method (look up “shooting methods”)

2 Likes

I’ve give the BVP problem a try but when the enery is not the eigen-energy, the wavefunction can be very hash. I don’t know if there’s any idea to deal with this like shooting and how to judge if the shape of wave function is resonable. Thanks for your quick comment!