SIRK schemes and Lokta-Volterra test

Question

I'm aiming at making a 2D compressible Euler equations numerical model, and to do so I need to implement a so-called "semi-implicit fourth order runge-kutta scheme" in order to integrate an equation of the form du/dt = f(u) + g(u); where g is a stiff term whereas f is not. Therefore, I've been following the scheme described by the following article https://www.researchgate.net/publication/228655062_Implicit-explicit_Runge-Kutta_method_for_combustion_simulation.

After quite some time spent on trying to understand how this works, I'm currently trying to test my newly implemented scheme with simple use cases. Following the example given on the scipy.integrate.solve_ivp documentation page, the Lokta-Volterra case, I'm trying to compare the results from my "SIRK45" scheme to the "RK45" and "Radau" schemes of solve_ivp (since my SIRK45 scheme requires two function, I take the other one as a null function).

While if the time integration step is small enough (~t_end/200) the explicit part of my scheme seems to behave fine. However for the same time step if we use the implicit part to do the integration, it looks like information is lost over time somehow.

The way I approached the implicit integration was to treat f(y + h*sum) as a constant, and rely on scipy.optimize.root function to solve iteratively equations of the type x - cst1 - g(y + dt*sum + dt*cst2*x) = 0. I believe there might be an issue in the way I treat the implicit part, but I can't figure out where or what.

---

Below is the code for the SIRK45 scheme I implemented (I can also add the test files codes should someone find them relevant).

import numpy as np
from scipy.optimize import root
def SIRK45(F,G,y0,h0,t0,tf,args,method='krylov'):
    #https://www.researchgate.net/publication/228655062_Implicit-explicit_Runge-Kutta_method_for_combustion_simulation
    ### Coefficients
    A = np.array([
        [1/4, 0, 0, 0, 0],
        [0.34114705729739, 1/4, 0, 0, 0],
        [0.80458720789763, -0.07095262154540, 1/4, 0, 0],
        [0.52932607329103, 1.15137638494253, -0.80248263237803, 1/4, 0],
        [0.11933093090075, 0.55125531344927, -0.1216872844994, 0.20110104014943, 1/4],
    ])
    B = np.array([0.11933093090075, 0.55125531344927, -0.1216872844994, 0.20110104014943, 1/4])
    E = np.array([
        [0, 0, 0, 0, 0],
        [0.39098372452428, 0, 0, 0, 0],
        [1.09436646160460, 0.33181504274704, 0, 0, 0],
        [0.14631668003312, 0.69488738277516, 0.46893381306619, 0, 0],
        [-1.33389883143642, 2.90509214801204, -1.06511748457024, 0.27210900509137, 0],
    ])
    
    ### Initialize
    y = y0
    t = t0
    h = h0
    
    ### Predictors
    print('Starting')
    f_k1 = F(y,*args)
    g_k1 = lambda x: x - f_k1 - G(y + h*(A[0,0]*x),*args)
    k1sol = root(g_k1, y, method=method)
    k1 = k1sol.x
    
    #
    print('Checkpoint 1/5')
    f_k2 = F(y + h*(E[1,0]*k1),*args1)
    g_k2 = lambda x: x - f_k2 - G(y + h*(A[1,0]*k1 + A[1,1]*x),*args)
    k2sol = root(g_k2, y, method=method)
    k2 = k2sol.x
    #
    print('Checkpoint 2/5')
    f_k3 = F(y + h*(E[2,0]*k1 + E[2,1]*k2),*args)
    g_k3 = lambda x: x - f_k3 - G(y + h*(A[2,0]*k1 + A[2,1]*k2 + A[2,2]*x),*args)
    k3sol = root(g_k3, y, method=method)
    k3 = k3sol.x
    #
    print('Checkpoint 3/5')
    f_k4 = F(y + h*(E[3,0]*k1 + E[3,1]*k2 + E[3,2]*k3),*args)
    g_k4 = lambda x: x - f_k4 - G(y + h*(A[3,0]*k1 + A[3,1]*k2 + A[3,2]*k3 + A[3,3]*x),*args)
    k4sol = root(g_k4, y,method=method)
    k4 = k4sol.x
    #
    print('Checkpoint 4/5')
    f_k5 = F(y + h*(E[4,0]*k1 + E[4,1]*k2 + E[4,2]*k3 + E[4,3]*k4),*args)            
    g_k5 = lambda x: x - f_k5 - G(y + h*(A[4,0]*k1 + A[4,1]*k2 + A[4,2]*k3 + A[4,3]*k4 + A[4,4]*x),*args)
    k5sol = root(g_k5, y, method=method)
    k5 = k5sol.x
    
    print('Finished, good luck pilot')
    ### Update the state vector
    new_y = y + h*(B[0]*k1 + B[1]*k2 + B[2]*k3 + B[3]*k4 + B[4]*k5)
    t = t+h0
    return t,new_y

EDIT - below the code to retrieve each SIRK45 plots :

def lotkavolterraCP(z, a, b, c, d):
    x, y = z
    return [a*x - b*x*y, -c*y + d*x*y]

### Useful
def zero_func(z, a, b, c, d):
    return np.array([0,0]).ravel()

### Numerical parameters
t = 0
tf = 15
nbr = 300
dt = 15/nbr
args=(1.5, 1, 3, 1)

y0 = np.array([10,5])
q1 = []
q2 = []

while t<tf:
    ### Call the semi-explicit scheme for integration
    # t,y = SIRK45(lotkavolterraCP,zero_func,y0.ravel(),dt,t,tf,args) #Explicit integration
    t,y = SIRK45(zero_func,lotkavolterraCP,y0.ravel(),dt,t,tf,args) #Implicit integration
    
    ### Update the parameters
    y0 = y
    q = y.reshape([2])
    q1.append(q[0])
    q2.append(q[1])

### Plotting
t = np.linspace(0, 15, nbr)
import matplotlib.pyplot as plt
plt.figure()
plt.plot(t, q1)
plt.plot(t, q2)
plt.xlabel('t')
plt.legend(['x', 'y'], shadow=True)
plt.title('Lotka-Volterra System - SIRK45')
plt.show()