必须实现一个系统ODE函数

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是方法、公差等的字典，请参阅文档。然后评估解决方案。

谢谢你Lutz Lehmann，我根据你的答案解决了这个问题。顺便说一句，方程3.3不是多余的，我应该澄清一下，它是一个表示每枚硬币中银含量的术语。所以它是一个独立的值（我同意论文中定义它的奇怪方式）。
以下是我最终使用的解决方案：

def dSdt(t,S):
x,y,z1,z2,z1_z2 = S

xn = p0*((z1_z2/p1)-1)

yn = p2*((x/p3)-1) + p4*((z1_z2/p1)-1) 

z1n = p7*z1*(1-(y/p6))

z2n = z2*(1-(y/p6))*(p7-p5*z1_z2)

z1_z2n = p5*(1-(y/p6))

return [xn,yn,z1n,z2n,z1_z2n]

# Before td

t1 = np.linspace(t0,395,abs(t0)+395,dtype=int)
S0_1 = [x_t0,y_t0,z1_t0,z2_t0,z1_div_z2_t0]
sol1 = odeint(dSdt,y0=S0_1,t=t1,tfirst=True)

# Handle division at td

sol1[-1][0] *= λ 
sol1[-1][1] *= λ
S0_2 = sol1.T[:,-1]

# After td

t2 = np.linspace(start=396,stop=500,num=500-396,dtype=int)
sol2 = odeint(dSdt,y0=S0_2,t=t2,tfirst=True)
sol= np.concatenate((sol1,sol2))
t = np.concatenate((t1,t2))