Как использовать жесткий решатель ODE в pyhton?

Question

Я новичок в программировании. Может ли кто-нибудь помочь мне решить эту жесткую систему ODE в python. Я проверил некоторые решения как для жесткой, так и для нежесткой системы. Но если я использую odeint в Python (так же, как ode45 в Matlab) и ode23s в Matlab, я получаю разные графики. Вот мой код в python с использованием odeint:

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

Vp = 3
Vi = 11
Vg = 10
E  = .2
tp = 6
ti = 100
td = 13
k  = 0.5
Rm = 209
a1 = 6.6
C1 = 300
C2 = 144
C3 = 100
C4 = 80
C5 = 26
Ub = 72
U0 = 4
Um = 94
Rg = 180
a  = 7.5
b  = 1.772
G = 0
kappa = (1/C4)*((1/Vi)-(1/(ti*E)))

TT=np.arange(1,51,1)
Y=len(TT)
d0=10**(-8)
A=5000
N=50
LL=np.zeros(Y)

def f1(x):
    return Rm/(1 + np.exp((-x/(Vg*C1))+a1))

def f2(x):
    return Ub*(1 - np.exp((-x/(C2*Vg))))

def f3(x):
    d1= (kappa*x)**(b)
    return (1/(C3*Vg))*((U0+(Um*d1))/(1+d1)) 

def f4(x):
    return Rg/(1+np.exp(a*((x/(C5*Vp))-1)))

def IG(x):
    return G

def ultradian(y,t):
    return f1(y[2])-E*((y[0]/Vp)-(y[1]/Vi))-(y[0]/tp),\
        E*((y[0]/Vp)-(y[1]/Vi))-(y[1]/ti),\
        f4(y[5])+IG(t)-f2(y[2])-f3(y[1])*y[2],\
        (1/td)*(y[0]-y[3]),\
        (1/td)*(y[3]-y[4]),\
        (1/td)*(y[4]-y[5])

for i in range(0,Y):
    T=TT[i]
    lya=np.zeros(N-1)
    y11_init = 0
    y12_init = 0
    y13_init = 0
    y14_init = 0
    y15_init = 0
    y16_init = 0

    y21_init = d0
    y22_init = 0
    y23_init = 0
    y24_init = 0
    y25_init = 0
    y26_init = 0


    for j in range(1,N):
        y1 = odeint(ultradian,[y11_init,y12_init,y13_init,y14_init,y15_init,y16_init],[0,T])
        y2 = odeint(ultradian,[y21_init,y22_init,y23_init,y24_init,y25_init,y26_init],[0,T])

        y1[-1][2]=y1[-1][2]+5000
        y2[-1][2]=y2[-1][2]+5000

        d=np.sqrt((y1[-1][0]-y2[-1][0])**2+(y1[-1][1]-y2[-1][1])**2+(y1[-1][2]-y2[-1][2])**2+
              (y1[-1][3]-y2[-1][3])**2+(y1[-1][4]-y2[-1][4])**2+(y1[-1][5]-y2[-1][5])**2)

       lya[j-1]=np.log(d/d0)
       y21new = y1[-1][0] + (d0/d)*(y2[-1][0]-y1[-1][0])
       y22new = y1[-1][1] + (d0/d)*(y2[-1][1]-y1[-1][1])
       y23new = y1[-1][2] + (d0/d)*(y2[-1][2]-y1[-1][2])
       y24new = y1[-1][3] + (d0/d)*(y2[-1][3]-y1[-1][3]) 
       y25new = y1[-1][4] + (d0/d)*(y2[-1][4]-y1[-1][4]) 
       y26new = y1[-1][5] + (d0/d)*(y2[-1][5]-y1[-1][5])

       y11_init = y1[-1][0]
       y12_init = y1[-1][1]
       y13_init = y1[-1][2]
       y14_init = y1[-1][3]
       y15_init = y1[-1][4]
       y16_init = y1[-1][5]

       y21_init = y21new
       y22_init = y22new
       y23_init = y23new
       y24_init = y24new
       y25_init = y25new
       y26_init = y26new

    LL[i]=np.mean(lya)

Вот изображение для Python, а также если я использую ode45 в Matlab: введите описание изображения здесь А вот изображение, если я использую ode23s в Matlab: введите описание изображения здесь