使用scipy.solve_ivp求解微分方程时结果随t_span起始值异常变化的问题求助

我在用solve_ivp求解微分方程，但发现结果会随问题设置变化出现奇怪的情况，我觉得这可能是求解器实现的问题而非代码问题，希望能得到大家的建议。

我的代码如下：

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

omega_ir = 1.0 * 2 * np.pi
gamma_ir = 0.5


def pulse(t):
    E0 = 1
    w = 0.35
    return -E0 * np.exp(-((t - 0) ** 2) / (w**2)) * np.sin((t - 0) / 0.33)


def equations(t, y):
    Q_ir, P_ir = y

    # Equations of motion
    dQ_ir = P_ir
    dP_ir = -(omega_ir**2) * Q_ir - gamma_ir * P_ir + pulse(t)

    return [dQ_ir, dP_ir]


initial_conditions = [0, 0]

# Time span for the simulation
t_span = (-1, 40)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

solution = solve_ivp(
    equations,
    t_span,
    initial_conditions,
    t_eval=t_eval,
)

Q_ir = solution.y[0]
print(solution.message)


fig, ax = plt.subplots()
ax.plot(t_eval, Q_ir / max(Q_ir))
ax.plot(t_eval, pulse(t_eval) / max(pulse(t_eval)))
ax.set_xlabel("Time")
ax.set_ylabel("Normalised intensity Intesnity")

plt.show()

这是一个简单的振子（类似单摆）问题：

• 当我运行上述代码时，一切正常，求解器返回的消息是：The solver successfully reached the end of the integration interval.
• 对应的结果图：蓝色曲线是归一化后的振子位置，橙色是归一化后的初始激励，曲线趋势合理——振子先跟随激励响应，之后呈现阻尼振荡逐渐衰减。

但当我把t_span修改为(-15, 40)时，结果就出现了异常：

• 此时蓝色曲线的实际数值量级约为1e-50（归一化后几乎贴合横轴），完全不符合预期。
• 我一开始以为是采样密度的问题，尝试把t_eval的采样点数增加到10000，但没有任何改善。而且只要t_span的起始时间早于-11，就会出现这个问题。

我猜测这可能是数学层面或者求解器实现的问题，也试过切换其他求解器方法，但得到的结果依然类似，都是近乎0的无意义值。

我对数值求解器的理论了解不多，希望有人能指点我解决这个问题的方向，或者告诉我这个问题是否无解，谢谢。

备注：内容来源于stack exchange，提问作者mmonti