[Python] SciPy solve_ivp cannot integrate over a derivative maximum in Friedmann's equations

I would like to numerically integrate the Friedmann equations to solve for the scalefactor in various universes:

import numpy as np
import scipy as sp

def H(t, a, om, ol, orad):
ok = 1 - om - ol - orad
return H0y * np.sqrt(np.maximum((orad / a**4) + (om / a**3) + (ok / a**2) + ol, 0.0))


def da_dt(t, a, om, ol, orad):
return a * H(t, a, om, ol, orad)

ol = -1.0
om = 1.0
orad = 0.02

t_span = [0, 100]
a0 = [1e-4]

sol = sp.integrate.solve_ivp(da_dt, t_span, a0, method='DOP853', args=(om, ol, orad), atol=1e-12, rtol=1e-10, max_step=0.05)

a_vals = sol.y[0]
t_vals = sol.t


This solver runs without throwing an error message (although I got errors with more lenient error tolerance and with some of the worse methods like RK45) but does not integrate correctly. For the parameters given, a plot of scalefactor (a) against time (t) should show an initial expansion followed by a contraction. Instead, it seems the solver gets stuck around da/dt = 0 and a(t) flatlines, staying the same value for all further values of t.

My friend tells me that MATLAB's ode45 can handle this ODE just fine and I'm failing to see why solve_ivp can't. If someone could provide insights, that would be much appreciated.

If someone wants the plotting code as well, here it is:

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 8))
plt.plot(t_vals, a_vals, label='$(\Omega_M,\Omega_\Lambda)$=(%.1f,%.1f)' % (om, ol))
plt.xlabel('Lookback Time (Gyr)')
plt.ylabel('Scalefactor')
plt.ylim(0, 3)
plt.legend(loc='upper left', frameon=False)
plt.show()

Continue reading...