Scipy solve_ivp extremely slow/freezing

I have been working on a problem that involves solving 2 pairs of coupled ode's. It has been recommended that we use scipy solve_ivp to do this however it seems to be either freezing or running extremely slowly (I have ran this for up to half an hour with no results).

I have limited experience with solving ode's numerically and have never used solve_ivp before so this could be a simple error but I have been unable to find it.

import numpy as np
from matplotlib import pyplot as plt
from scipy import integrate as sc

param_dict = {
    'r1': 0.36,
    'r2': 0.15,
    'k1': 4000000,
    'b12': 0.0000002,
    'b21': 0.69,
    'p1': 80,
    'p2': 275,
    'q1': 0.0002,
    'q2': 0.0004,
    'phi1': 0.00095,
    'phi2': 0.000075,
    'c1': 20000,
    'c2': 30000
}

# INITIAL VALUES
init_X = np.array([4000000,275000])
init_E = np.array([800,200])
y0 = np.concatenate((init_X, init_E))

def functions(t, X, r1, r2, k1, b12, b21, q1, q2, phi1, phi2, c1, c2, p1, p2):
    x1, x2, y1, y2 = X
    dx1dt = (r1 * x1) * (1 - (x1 / k1) - (b12 * x2)) - (q1 * y1 * x1)
    dx2dt = (r2 * x2) * (1 - (x2 / (b21 * x1))) - (q2 * y2 * x2)
    dy1dt = phi1 * y1 * (p1*q1 * x1 - c1)
    dy2dt = phi2 * y2 * (p2*q2 * x2 - c2)
    return [dx1dt, dx2dt, dy1dt, dy2dt]

t_span = [0,20]
t = np.linspace(0,100,1000)

res = sc.solve_ivp(functions, t_span, y0, args=tuple(param_dict.values()), dense_output=True)

zz = res.sol(t)
x1 = []
x2 = []
e1 = []
e2 = []
for ii in range(len(zz.T)):
    x1.append(zz.T[ii,0])
    x2.append(zz.T[ii,1])
    e1.append(zz.T[ii,2])
    e2.append(zz.T[ii,3])

plt.plot(x1)
plt.plot(x2)

If I set init_E = [0,0] it runs seemingly as expected however any other value here seems to lead to this freezing/extremely slow problem.

The model here describes a predator/prey model for fisheries management for those interested.

Solution

I think the problem might be that the keys in your dictionary are not in the same order as the arguments to your function.

If I do the following it solves immediately.

param_names = ['r1', 'r2', 'k1', 'b12', 'b21', 'q1', 'q2', 'phi1', 
               'phi2', 'c1', 'c2', 'p1', 'p2']
args = tuple(param_dict[p] for p in param_names)
res = sc.solve_ivp(functions, t_span, y0, args=args, dense_output=True)

(Not sure if this is the correct solution but it doesn't hang up).

Bit of advice:

In general it's always good to make at least one call to your function and verify its output before passing it to a solver. E.g. at the initial condition like this:

f0 = functions(0, y0, *args)
assert np.array_equal(f0, [-719200.0, 15139.945652173912, 33440.0, 3.75])