Setting `abstol` and `reltol` when solving the Schrödinger equation with OrdinaryDiffEq

[Mason]

I was curious and tried this out using Magnus solvers and fixed time grid as you mentioned you needed that, and it seems to have performed very well:

using OrdinaryDiffEq, SciMLOperators, Plots

function tst2(;)
    tlist = range(-100, 100; length=1001)
    
    Ψ₀ = [-0.999998+0.0im, 0.0021739+0.0im]
    
    function update_func!(W, u, p, t)
        Δ = 1.0
	    α = 4.6
        W[1,1] = -im*(0.5 * α * t)
        W[2,1] = -im*(Δ)
        W[1,2] = -im*(Δ)
        W[2,2] = +im*(0.5 * α * t)
    end
    W = MatrixOperator(zeros(ComplexF64, 2, 2); update_func!) # -im * H
    
    prob = ODEProblem(W, Ψ₀, (first(tlist), last(tlist)))
    
    Ψ = solve(prob, MagnusMidpoint(), dt=step(tlist), saveat=tlist)

    popdyn_ode_direct = reduce(hcat, map(t -> abs2.(Ψ(t)), tlist))
    p = plot(tlist, popdyn_ode_direct'; label=["|⟨Ψ(t)|ϕ₀⟩|²" "|⟨Ψ(t)|ϕ₁⟩|²"], legend=:right, title="Projection onto canonical basis states", xlabel="time", ylabel="population")
    plot!(tlist, norm.(Ψ.(tlist)); label="|⟨Ψ(t)|Ψ(t)⟩|")
    p
end

Here’s a zoom-in on the drift of |⟨Ψ(t)|Ψ(t)⟩| and it is very precise:

This solver is also extremely fast so that’s nice (8x faster than Tsit5() for this example)