solve_ivp gets stuck and I can't figure out why

I am using SciPy's solve_ivp to simulate a model consisting of 5 equations:

Upsilon = 0.305
Alpha = 0.00000055
r = 0.176
b = 0.0000000000005

def model(t, y0, Beta, Eta, Rmin, p1, p2, p3, A, a, Xi, Epsilon, Lambda, Theta, Mu, Delta, Gamma):

    CD, CT, CM, CE, T, TOTAL = y0

    def k_t(t):
        return Rmin + (p1 / (1 + (p2 * t) ** p3))

    def F_t(T):
        return T / (A + T)

    def f_t(CF, T):
        return (CF / T) / (Upsilon + ((a + CF) / T))

    dCD_dt = -(Beta + Eta) * CD
    dCT_dt = (Eta * CD) + (k_t(t) * F_t(T) * CT) - ((Xi + Epsilon + Lambda) * CT) + (Theta * T * CM) - (Alpha * T * CT)
    dCM_dt = (Epsilon * CT) - (Theta * T * CM) - (Mu * CM)
    dCE_dt = (Lambda * CT) - (Delta * CE)
    dT_dt = ((r * T) * (1 - (b * T))) - (Gamma * f_t((CD + CT), T) * T)

    TOTAL = dCD_dt + dCT_dt + dCM_dt + dCE_dt

    print(t)
    return [dCD_dt, dCT_dt, dCM_dt, dCE_dt, dT_dt, TOTAL]

I am performing sensitivity analysis, so I am running the model about 37,000 times with variations in the input parameters. When I run the model, I do so with the following code:

#p_values is a numpy array

final = 100
t_range = (0, final)
t_eval = np.arange(0, final, 0.1)

# run model for each parameter combination
for i in range(len(p_values)):
    CD = p_values[i][15]
    CT, CM, CE = 0, 0, 0
    T = 10 ** 7
    TOTAL = CD
    y0 = [CD, CT, CM, CE, T, TOTAL]

    solution = solve_ivp(model, t_range, y0, t_eval=t_eval, args=(p_values[i][0], p_values[i][1], p_values[i][2], p_values[i][3], p_values[i][4], p_values[i][5], p_values[i][6], p_values[i][7], p_values[i][8], p_values[i][9], p_values[i][10], p_values[i][11], p_values[i][12], p_values[i][13], p_values[i][14])).y

But at some point while it's running it will just... get stuck. Sometimes it happens early (within the first thousand solutions) and sometimes it can get to 10,000 solutions or more before it gets stuck.

To try and solve the problem, I started printing t for each loop of the model, and I noticed something very strange - when it gets stuck, it seems that t stays stagnant, or at least grows extremely slowly. This doesn't make sense to me, as I assume it should always increase by 0.1 at each step. Whenever this happens t is quite high (above 95), but I don't know if this is a coincidence or not. If not, I suppose it could be an issue with the function k(t) in the model producing a number that is far too small, but I don't see why this should be an issue.

If anyone with more experience than me has a suggestion to try and fix it, that would obviously be much appreciated. In lieu of that, I would also be okay with having some way to skip over that run of the model and move onto the next when the running time gets too long, but I don't know how to implement something like this without making the runs very slow.

You could possibly try using the optional method parameter in your solve_ivp, and set it to method = "BDF". It might have a little more accuracy.

According to
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html :

> ‘BDF’: Implicit multi-step variable-order (1 to 5) method based on a
> backward differentiation formula for the derivative approximation [5].
> The implementation follows the one described in [6]. A quasi-constant
> step scheme is used and accuracy is enhanced using the NDF
> modification. Can be applied in the complex domain.