The most immediate fix is to lower the relative tolerance:
ode_sol_direct = ode_solve(ode_problem_direct, ODE.Tsit5(); abstol=1e-12, reltol = 1e-10, saveat=tlist)
Typically, although I don’t think it affects the accuracy, there’s also no need to include the saveat = tlist. You can get get the values of the solution by interpolation by calling ode_sol_direct(t) after you’ve solved it.
This is not my area, but Tsit5() is probably not the best choice for solver for this problem, check the “Probelm Types” heading in the left bar of DifferentialEquations.jl: Efficient Differential Equation Solving in Julia · DifferentialEquations.jl . I don’t remember Quantum Mechanics stuff well enough to figure out if any of those types would match up with your system.