Issue solving differential equation with solve_ivp depending on t_span

I'm using solve_ivp but I'm getting strange results depending on the settings of the problem, though i think it's more an issue with the implementation of the solver than a coding issue, but I'd love for someone to provide an input. So my code is as follows:

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()

So it's a simple pendulum/oscillator problem. if I run this it works without issues, the message is: The solver successfully reached the end of the integration interval. great. The plot is as follows (blue is the position of the oscillator, orange is the initial excitation, both normalised for clarity): Very nice. However if I try to change the t_span=(-15, 40) The plot is as follows: Note that is still normalised, the blue line is actually ~1e-50 or something like that. I of course thought it was an issue with sampling, so I tried to change it to finer sampling (like 10000 points) without success. This seem to happen if the first time point is earlier than -11.

I think this is an issue with the math or the implementation of the solver. I tried changing to other methods, but it seems to give similar nonsensical (~0) results. I don't know much about the theory of these solvers, so I'd be happy if someone could point me to the direction of solving this issue, or let me know if this is not solvable. Thanks

Solve_ivp uses an adaptive time step. This is based on how fast things are happening. If you start from well before any pulse then nothing appears to happen, the solver rapidly increases the time step and then … misses the pulse when it does occur!

A simple solution is to limit the time step in solve_ivp. Add the optional parameter max_step=0.1 (or anything else sufficiently small) and you should be fine.