Solving ODEs python with non-independent funcitons python

Question

I'm trying to plot a multi-equation ODE from this paper and I'm having trouble figuring out how to express the dependencies between the functions. I.e, in 3.1, it uses both z1 and z2. Aswell as using its own precious state x. I've checked the Scipy documentation but they don't seem to express the functionality that I'm trying to achieve.

How do I express I.e equation 3.1 with Scipy?

Answer 1

You have to implement a system ODE function

def sys(t,u):
    x,y,z1,z2 = u
    v1 = z1/(p1*z2)-1
    v2 = 1-y/p6
    return p0*v1, p2*(x/p3-1)+p4*v1, p7*z1*v2, (p7-p5*z2/z1)*z2*v2

#integrate to jump
sol1 = solve_ivp(sys,(t0,td),u0,**options)
# implement the jump
x,y,z1,z2 = sol1.y[:,-1]
fac = (2-λ)/(2+λ); 
x*=fac; y*=fac
ud = [x,y,z1,z2]
# restart integration
sol2 = solve_ivp(sys,(td,tf),ud,**options)

options is a dictionary for the method, tolerances, etc., see documentation. Then evaluate the solutions.