是什么导致PythonOdeint在某一点之后返回奇怪的结果？

Question

我对python odeint函数有问题。我正在使用扩展的SIR模型来模拟传染病的传播，并使用python来求解耦合微分方程，然后绘制它们。通常情况下，这些图看起来非常正常，并返回t的第一个x值的合理和准确的结果，但之后它们膨胀到无穷大，或者创建明显不正确的奇怪模式。我对python和odeint都比较陌生，所以我很难找到问题所在。这是我的密码：

def plotsirhd(t, S, I, R, H, D):
    f, ax = plt.subplots(1,1,figsize=(10,4))
    ax.plot(t, S, 'cornflowerblue', alpha=0.7, linewidth=2, label='Susceptible')
    ax.plot(t, I, 'r', alpha=0.7, linewidth=2, label='Infected')
    ax.plot(t, R, 'Maroon', alpha=0.7, linewidth=2, label='Hospitalised')
    ax.plot(t, H, 'darkslateblue', alpha=0.7, linewidth=2, label='Recovered')
    ax.plot(t, D, 'k', alpha=0.7, linewidth=2, label='Dead')
    ax.set_xlabel('Time (days)')
    ax.set_ylabel('Fraction of population',labelpad=14)

    ax.yaxis.set_tick_params(length=0)
    ax.xaxis.set_tick_params(length=0)
    ax.grid(b=True, which='major', c='w', lw=2, ls='-')
    legend = ax.legend(borderpad=2.0)
    legend.get_frame().set_alpha(0.5)
    for spine in ('top', 'right', 'bottom', 'left'):
        ax.spines[spine].set_visible(False)
    plt.show();

    plt.show();

def diffeq3(y, t, beta3, q3, gamma=1/9, delta=1/7, aitch=1/12, p=0.05, q0=0.45, q1=1, c=2/6665, L0=46, L1=120, L=0.8):
    def q3(t, q0, q1, c):
        if H < c:
            return q0
        else:
            return ((H-c)*q1+c*q0)/H
    def beta3(t, L0, L):
        if t < L0:
            return gamma * (1.8/(1+np.exp(-0.5*(-t+L0)))+1.8)
        else:
            if L0<=t<=L1:
                return L * gamma
            else:
                if H > 3*c/4:
                    return L * gamma
                else:
                    return 1.8 * gamma
    S, I, R, H, D = y
    dSbydt = -1*beta3(t, L0, L)*S*I
    dIbydt = beta3(t, L0, L)*S*I - (p*aitch*I + gamma*(1-p)*I)
    dHbydt = p*aitch*I - (q3(t, q0, q1, c)*delta*H + (1-q3(t, q0, q1, c))*delta*H)
    dRbydt = (1-p)*gamma*I + (1-q3(t, q0, q1, c))*delta*H
    dDbydt = q3(t, q0, q1, c)*delta*H
    return dSbydt, dIbydt, dRbydt, dHbydt, dDbydt


t = np.linspace(0, 899, 900) 
y0 = S0, I0, H0, R0, D0 
gamma=1/9
delta=1/7
aitch=1/12
p=0.05
q0=0.45
q1=1
c=2/6665
L0=46
L1=120
L=0.8
ret = odeint(diffeq3, y0, t, args=(beta, q, gamma, delta, aitch, p, q0, q1, c, L0, L1, L))
S, I, H, R, D = ret.T

plotsirhd(t, S, I, R, H, D)

此外，如果经常多次重新运行代码，则会出现不同的错误。对于我的beta函数和不同参数的细微变化，这个问题仍然存在。有没有办法解决这个问题？我希望看到耦合函数的输出，以获得更大的t值。多谢各位