How to use a spike firing event in scipy.integrate.solve_ivp?

Question

I have tried to simulate an adaptive exponential integrate-and-fire neuron using the scipy.integrate.solve_ivp method. In the model there is a resetting condition which I was successfully able to implement by defining an event function which resets the value of membrane potential when it crosses a particular threshold. It works well except for the cases when the reset potential is greater than the the threshold potential. Instead of registering a spike, it continues to grow according to the equation beyond the event condition.

This is the code I have written for the simulation:

# Importing packages
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, odeint
import pandas as pd

# Parameters
adaptation = False  # When the reset is below threshold (works fine)
regular_bursting = True  # When the reset is above threshold

if adaptation:
    C = 200.0  # Membrane capacitance (pF)
    g_l = 12.0  # Leak conductance (nS)
    E_l = -70.0  # Leak reversal potential (mV)
    V_t = -50.0  # Threshold potential (mV)
    V_peak = 20.0  # Peak potential (mV)
    delta_t = 2  # slope factor (mV)
    tau_w = 300.0  # Adaptation time constant (ms)
    a = 2.0  # Subthreshold adaptation coupling (nS)
    b = 60  # Spike-triggered adaptation increment (pA)
    V_reset = -58.0  # Reset current
    current = 500  # pA

if regular_bursting:
    C = 200.0  # Membrane capacitance (pF)
    g_l = 10.0  # Leak conductance (nS)
    E_l = -58.0  # Leak reversal potential (mV)
    V_t = -50.0  # Threshold potential (mV)
    V_peak = 20.0  # Peak potential (mV)
    delta_t = 2  # slope factor (mV)
    tau_w = 120.0  # Adaptation time constant (ms)
    a = 2.0  # Subthreshold adaptation coupling (nS)
    b = 100  # Spike-triggered adaptation increment (pA)
    V_reset = -46  # Reset current
    current = 210  # pA


# Simulation conditions
t_sim = 1000
dt = 0.05


# Defining spiking event
def spike_event(t, state):

    return state[0] - V_t  # Returns true if event is hit

spike_event.terminal = True  # Stops integration if event is hit


# Initial conditions
state_ini = [E_l, 0]


# Differential equation
def adex_neuron(t, state):
    V, w = state

    I = current
    dV_dt = (g_l * (E_l - V) + g_l * delta_t * np.exp((V - V_t) / delta_t) - w + I) / C 
    dw_dt = (a * (V - E_l) - w) / tau_w  

    return [dV_dt, dw_dt]


# Solving ode
time_dump = []
potential_dump = []
t = 0 # initial time

spike = [] # Stores spike times

while t < t_sim:
    sol = solve_ivp(adex_neuron, (t, t_sim), y0=state_ini, events=spike_event, method="RK45")
    time_dump.append(sol.t)    
    potential_dump.append(sol.y[0])

    if sol.status == 1:  # If event is hit
        t = sol.t[-1]
        spike.append(sol.t_events[0][0])
        state = sol.y[:, -1].copy()
        state[0] = V_reset
        state[1] += b
        state_ini = state
    else:
        break


plt.plot(np.concatenate(time_dump)[:], np.concatenate(potential_dump)[:])
plt.xlabel("Time")
plt.ylabel("Potential")
plt.show()

The V vs time graph - It can be seen the graph goes beyond the threshold value:

The initial condition after the reset value is not being considered by the event function. Any help would be appreciated.