It turns out that the code within evolution_rate, i.e. the Numpy backend, was correct, but yielded unphysical values and a FloatingPointError when dt=time_step was 1/number_of_steps = 1/1e4 = 1e-4.
So the error occurred because the time step was too large, not because the computation was incorrect.
The Numba backend, which calculates numerically the same time derivatives, somehow surpasses this FloatingPointError.
Solved by setting number_of_steps = 1e6 instead of number_of_steps = 1e4, which will give a runtime of 2 hours for the Numpy backend and a runtime of 2 minutes for the Numba backend.
It turns out that the code within
evolution_rate, i.e. the Numpy backend, was correct, but yielded unphysical values and a FloatingPointError whendt=time_stepwas1/number_of_steps = 1/1e4 = 1e-4. So the error occurred because the time step was too large, not because the computation was incorrect. The Numba backend, which calculates numerically the same time derivatives, somehow surpasses this FloatingPointError. Solved by settingnumber_of_steps = 1e6instead ofnumber_of_steps = 1e4, which will give a runtime of 2 hours for the Numpy backend and a runtime of 2 minutes for the Numba backend.