[Python] graphing nonlinear systems of ODEs in python

the model Im trying to graph this model however I am just not sure weather my code is leading me to a reasonable solution:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp


K = 17.3
epsilon = 4.9
r = 2.7
mu = 0.9
gamma = 3.2
alpha = 4 / (K - epsilon) ** 2

def mosquito_odes(t, Y, alpha, r, K, epsilon, mu, gamma):
N, M = Y

dNdt = alpha * r * (N**2 / (N + M + 1e-6)) * (K - N) * (N - epsilon) - mu * N
dMdt = gamma * N / (1 + N) - mu * M

return [dNdt, dMdt]

N0 = K/1.5
M0 = 5.6667
Y0 = [N0, M0]

T = 5
t_eval = np.linspace(0, T, 100)


sol_ode = solve_ivp(
mosquito_odes,
[0, T],
Y0,
args=(alpha, r, K, epsilon, mu, gamma),
method="RK45",
t_eval=t_eval
)


t_vals = sol_ode.t
N_vals, M_vals = sol_ode.y

plt.figure(figsize=(8, 5))
plt.plot(t_vals, N_vals / K, label="Wild Population (N)", linestyle="-", color="b")
plt.plot(t_vals, M_vals / K, label="Sterile Population (M)", linestyle="--", color="r")
plt.axhline(1.9 / K, linestyle="dashed", color="brown", label="Threshold γc=1.9") # Show threshold
plt.xlabel("Time")
plt.ylabel("Fraction of Population")
plt.title("Non-Spatial Mosquito Population Dynamics")
plt.legend()
plt.grid()
plt.show()



the outcome is not what I expect (for these parameters I expect both ODEs to go to 0. however it seems that only one is.

Continue reading...