Inconsistency in ODE solution using solve_ivp when adding events

I am using scipy.solve_ivp to solve an ODE in a loop, such that in every iteration I change some term in the ODE (which is like the force applied).

In terms of my problem - I just cant seem to understand why when I track events of zero crossing the solution sometimes stalls and essentially zeros out (as the figure below shows). The code is a bit long, but only a small portion of it is actually relevant to my problem

(1) initialize some constants

MASS = 2.6e-4
WING_LENGTH = 0.07  # meters
GYRATION_RADIUS = 0.6 * WING_LENGTH
MoI = MASS * GYRATION_RADIUS ** 2 # Moment of Inertia
AIR_DENSITY = 1.2
WING_AREA = 0.5 * WING_LENGTH * (0.5 * WING_LENGTH) * np.pi  # 1/2 ellipse with minor radios ~ 1/2 major = length/2
# drag coefficients from whitney & wood (JFM 2010):
C_D_MAX = 3.4
C_D_0 = 0.4
C_L_MAX = 1.8
ZERO_CROSSING = 1
RADIAN45 = np.pi / 4
RADIAN135 = 3 * RADIAN45

(2) Define a wrapper which is a RobotSimulation class - it's an easy interface that allows to easily change the values of the ODE. The important function is solve_dynamics (the others are not part of the problem I believe)

def phi_dot_zero_crossing_event(t, y):
    """
    this event is given to solve_ivp to track if phi_dot == 0
    :param t: unused time variable
    :param y: a vector of [phi,phi_dot]
    :return:
    """
    return y[1]


class RobotSimulation:
    def __init__(self, motor_torque=lambda x: 0, alpha=RADIAN45, phi0=0.0, phi_dot0=0.01, start_t=0, end_t=0.05, delta_t=0.001):
        self.solution = 0
        self.motor_torque = motor_torque
        self.alpha = alpha  # angle of attack
        self.phi0 = phi0
        self.phi_dot0 = phi_dot0
        self.start_t = start_t
        self.end_t = end_t
        self.delta_t = delta_t

    def set_time(self, new_start, new_end) -> None:
        """
        we update the start time to be the end of the prev time and the end time
        as the last end time + the window size we wish
        """
        self.start_t = new_start
        self.end_t = new_end

    def set_init_cond(self, new_phi0, new_phi_dot0) -> None:
        """
        we update the initial conditions to be the last value of previous iterations
        """
        self.phi0 = new_phi0
        self.phi_dot0 = new_phi_dot0

    def set_motor_torque(self, new_motor_torque) -> None:
        """
        sets the moment to the new value
        """
        self.motor_torque = new_motor_torque

    def flip_alpha(self) -> None:
        """
        (Irrelevant for my question) changes the angle of attack's sign
        :return:
        """
        self.alpha = RADIAN135 if self.alpha == RADIAN45 else RADIAN45

    def c_drag(self) -> float:
        """
        calculates the drag coefficient based on the angle of attack
        """
        return (C_D_MAX + C_D_0) / 2 - (C_D_MAX - C_D_0) / 2 * np.cos(2 * self.alpha)

    def drag_torque(self, phi_dot):
        """
        the drag moment
        """
        return 0.5 * AIR_DENSITY * WING_AREA * self.c_drag() * GYRATION_RADIUS * phi_dot * np.abs(phi_dot)

    def phi_derivatives(self, t, y):
        """
        A function that defines the ODE that is to be solved: I * phi_ddot = tau_z - tau_drag.
        We think of y as a vector y = [phi,phi_dot]. the ode solves dy/dt = f(y,t)
        :return:
        """
        phi, phi_dot = y[0], y[1]
        dy_dt = [phi_dot, (self.motor_torque(t) - self.drag_torque(phi_dot)) / MoI]
        return dy_dt

    def solve_dynamics(self, phiarr, phidotarr, phiddotarr, angarr, torquearr):
        """
        solves the ODE
        :return:
        """
        phi_0, phi_dot_0 = self.phi0, self.phi_dot0
        start_t, end_t, delta_t = self.start_t, self.end_t, self.delta_t
        phi_dot_zero_crossing_event.terminal = True

        ang = []
        times_between_zero_cross = []
        sol_between_zero_cross = []
        while start_t < end_t:
            torquearr.append(self.motor_torque(start_t))
            sol = solve_ivp(self.phi_derivatives, t_span=(start_t, end_t), y0=[phi_0,phi_dot_0],events=phi_dot_zero_crossing_event)
            self.solution = sol.y
            ang.append(self.alpha * np.ones(len(sol.t)))  # set alpha for every t based on solution's size
            times_between_zero_cross.append(sol.t)
            sol_between_zero_cross.append(sol.y)
            if sol.status == ZERO_CROSSING:
                start_t = sol.t[-1] + delta_t
                phi_0, phi_dot_0 = sol.y[0][-1], sol.y[1][-1]  # last step is now initial value
                self.flip_alpha()
            else:  # no zero crossing
                break
        time = np.concatenate(times_between_zero_cross)
        phi, phi_dot = np.concatenate(sol_between_zero_cross, axis=1)
        ang = np.concatenate(ang)

        _, phi_ddot = self.phi_derivatives(time, [phi, phi_dot])
        phiarr.append(phi)
        phidotarr.append(phi_dot)
        phiddotarr.append(phi_ddot)
        angarr.append(ang)

(3) Finally the main function calls solve_dynamics sequentially for different values of actions and updates the time window and initial conditions:

phi_arr, phi_dot_arr, phi_ddot_arr, angarr, torquearr = [], [], [], [], []
phi0, phidot0 = 0, 0.01
start_t, end_t, delta_t = 0, 0.05, 0.001
sim = RobotSimulation()
for action in [1, -1, 1, -1, 1, -1]:
    sim.set_motor_torque(lambda x: action)
    sim.solve_dynamics(phi_arr, phi_dot_arr, phi_ddot_arr, angarr, torquearr)
    phi0, phidot0 = sim.solution[0][-1], sim.solution[1][-1]
    start_t = end_t
    end_t += 0.05
    sim.set_init_cond(phi0, phidot0)
    sim.set_time(start_t, end_t)
phi_arr = np.concatenate(phi_arr)
phi_dot_arr = np.concatenate(phi_dot_arr)
phi_ddot_arr = np.concatenate(phi_ddot_arr)
angle_arr = np.concatenate(angarr)

Intuitively, as the actions I apply are just repeating +1 and -1, I cannot understand why the blue points in the figure (see below) act so different