Problems with stability of Scipy.integrate solve_ivp

I'm using SciPy to solve a system of differential equations, and I am having some problems with the stability of the solution. Essentially, there are two state variables, and they should go to two relatively stable solutions, one asymptotically approaching zero (essentially zero) and the other to 0.59.

When I do this in R with the deSolve package it works fine, but I'm trying to use scipy and I'm getting a weird thing where the final solution, the state variables are not stable even though they should be? The values for the state variable that should be at zero go from a series of values in the order of 7.41e-323 and then jump up to 5.3e-001 for one step and then go back down. I have no idea why this is, but I'm wondering if anyone can provide any advice on how to a) fix this, or b) another option to use other than scipy?

Trying this in both R (lsoda package) and in Maple have yielded stable solutions.

Thanks!

Code:

# %%
import numpy as np
import scipy.integrate as spi
from scipy.integrate import solve_ivp

# %%
x_coords = np.arange(0.01,1,0.05)
#this makes a mesh grid at those points (the whole square)
M_temp, C_temp = np.meshgrid(x_coords, x_coords)

#only includes the points I want (the lower triangular in M x C space)
C_points = C_temp[M_temp + C_temp <=1]
M_points = M_temp[M_temp + C_temp <=1]
T_points = 1 - C_points - M_points
# initialize array
M_array = np.zeros((50000,len(C_points)))
C_array = np.zeros((50000,len(C_points)))

# %%
for i in range(0,len(C_points)):
    def rhs(s,v):
        a = 0.25
        g = 0.3
        z = 0.05
        r = 0.55 
        d = 0.24 
        y = 0.77 
        return [r*v[0]*v[2] + z*r*v[2] - d*v[0] - a*v[0]*v[1], a*v[0]*v[1] -
        (g*v[1])/(v[1]+v[2]) + y*v[1]*v[2],-r*v[0]*v[2] - z*r*v[2] + d*v[0]+
        (g*v[1])/(v[1]+v[2]) - y*v[1]*v[2]]

    res = solve_ivp(rhs, (0, 50000),
        [C_points[i], M_points[i], T_points[i]],
        t_eval =np.arange(0,50000,1)) #solves from t =0 -> t = 50000
    M_array[:,i] = res.y.T[:,1]
    C_array[:,i] = res.y.T[:,0] #[0:50,2]

trajectories = np.ones(len(C_points)*len(M_array[:,i]),\
    dtype = {'names': ['M', 'C'], 'formats':[np.float64, np.float64]})
trajectories['M'] = M_array.flatten()
trajectories['C'] = C_array.flatten()

print(trajectories['M'][1049900:1049999])

Output:

Solution

If you print the first decimal after the dot of the last M component in a C-M grid, you get

0.01: 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.06: 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.11: 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.16: 0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.21: 0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5  
0.26: 0 0 0 0 5 5 5 5 5 5 5 5 5 5 5  
0.31: 0 0 0 0 0 5 5 5 5 5 5 5 5 5  
0.36: 0 0 0 0 0 0 5 5 5 5 5 5 5  
0.41: 0 0 0 0 0 0 0 5 5 5 5 5  
0.46: 0 0 0 0 0 0 0 0 5 5 5  
0.51: 0 0 0 0 0 0 0 0 0 5  
0.56: 0 0 0 0 0 0 0 0 0  
0.61: 0 0 0 0 0 0 0 0  
0.66: 0 0 0 0 0 0 0  
0.71: 0 0 0 0 0 0  
0.76: 0 0 0 0 0  
0.81: 0 0 0 0  
0.86: 0 0 0  
0.91: 0 0  
0.96: 0

This pattern strongly suggests that the observed pattern is not a random error, but a property of the system, a certain set of initial values, roughly for C0 < M0, leads to a fixed point with a non-zero M component, 0.53371702 in the best estimate.

A streamplot for the first two components confirms this hypothesis, the saddle point at [0.31024394, 0.22563217] separates the regions converging to the stable points [0.07006604, 0.53371702] and [0.59734069, 0.0]

Check again that the code in both R and python versions really does the same. Find out the integration method in R, it is probably lsoda. Then find out the default tolerances, or set sensible tolerances (atol = 1e-11, rtol=1e-13) that then are replicated in python. Use the method(="BDF") option to set the integration method of solve_ivp.