For the same problem why NumbaLSODA gives Segmentation Fault but solve_ivp works fine?

I am solving a system of Odes using python. My attempt using solve_ivp method works fine. However, the same problem using NumbaLSODA gives segmentation fault. Can anybody please help?

################## solve_ivp code: ##################

import os,sys
import numpy as np
from numpy.linalg import multi_dot
from scipy import integrate,interpolate
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from tqdm import tqdm
import time
from scipy.interpolate import RegularGridInterpolator
import tarfile

mphi = 1.
GammaSM = 0.1*mphi 
gSM = 106.75
Mpl = 1. 


def fEQ(k,Y,x):
    fEQdenom = np.exp(np.sqrt(k**2 + (x*mphi)**2)/Y) - 1
    fEQval = pow(2*np.pi,-3)*(1./fEQdenom)
    return fEQval 

xipow = np.arange(-2.,3.+1,1)
xiarray = 10**xipow
increyxi = 1.

def system(t,u, m0, pbar, state):
    last_t, dt = state
    
    time.sleep(0.1)
    n = int((t - last_t)/dt)
    pbar.update(n)
    
    state[0] = last_t + dt * n    
    
    du = np.zeros((len(xipow)+1,))
    p = np.zeros((len(xipow),))

    rhophi = 0
    for ii in range(len(xipow)):
        rhophi += pow(xiarray[ii],3)*np.sqrt(pow(xiarray[ii],2)+(mphi*t)**2)*u[ii]  

    Hubsq = (1./(3*pow(Mpl*t*t,2)))*(4*np.pi*increyxi*rhophi + pow(np.pi,2)*(gSM/30.)*pow(u[len(xipow)],4))

    Hub = np.sqrt(Hubsq)

    for ii in range(len(xipow)):
        coeff = (mphi*GammaSM)/(Hub*np.sqrt(pow(xiarray[ii],2)+pow(t*mphi,2)))
        p[ii] =  fEQ(xiarray[ii],u[len(xipow)],t) - u[ii] 
        du[ii] = coeff *p[ii] 

    coefftemp = - ((30*mphi*GammaSM)/(np.pi*gSM*Hub))*pow(u[len(xipow)],-3)
    temprhs = np.sum(p*xiarray*xiarray*xiarray)
    du[len(xipow)] = coefftemp *temprhs * increyxi     
        
    return du


u0 = np.array([2.49282627e-12, 2.49282621e-12, 2.49282004e-12, 2.49220314e-12, 2.43131611e-12, 2.35213266e-13, 100])

ts1 = 1.
ts2 = 10.
ts = np.array([ts1,ts2])
times_step = np.linspace(ts[0],ts[1],90)

with tqdm(total=100, unit="‰") as pbar:
    
    sol = solve_ivp(system, ts, u0, method='LSODA', t_eval=times_step, rtol=1e-6, atol=1e-6, args=[1, pbar, [ts[0], (ts[1]-ts[0])/101]])

############# The same system of equations using NumbaLSODA: #############

from NumbaLSODA import lsoda_sig, lsoda
from numba import njit, cfunc
import numpy as np
import timeit
import numba as nb


mphi = 1. 
GammaSM = 0.1*mphi
gSM = 106.75
Mpl = 1. 

@njit
def fEQ(k,Y,x):
    fEQdenom = np.exp(np.sqrt(k**2 + (x*mphi)**2)/Y) - 1
    fEQval = pow(2*np.pi,-3)*(1./fEQdenom)
    return fEQval 
    
xipow = np.arange(-2.,3.+1,1)
xiarray = 10**xipow
increyxi = 1.


@cfunc(lsoda_sig)
def rhs(t, u, du, p):
    u_ = nb.carray(u, (len(xipow)+1,))
    p_ = nb.carray(p, (len(xipow),))

    rhophi = 0
    for ii in range(len(xipow)):
        rhophi += pow(xiarray[ii],3)*np.sqrt(pow(xiarray[ii],2)+(mphi*t)**2)*u[ii]  

    Hubsq = (1./(3*pow(Mpl*t*t,2)))*(4*np.pi*increyxi*rhophi + pow(np.pi,2)*(gSM/30.)*pow(u[len(xipow)],4))

    Hub = np.sqrt(Hubsq)

    temprhs = 0

    for ii in range(len(xipow)):
        coeff = (mphi*GammaSM)/(Hub*np.sqrt(pow(xiarray[ii],2)+pow(t*mphi,2)))
        p[ii] =  fEQ(xiarray[ii],u[len(xipow)],t) - u[ii] 
        temprhs += pow(xiarray[ii],3)*p[ii] 
        du[ii] = coeff * p[ii]

    coefftemp = - ((30*mphi*GammaSM)/(np.pi*gSM*Hub))*pow(u[len(xipow)],-3)
    du[len(xipow)] = coefftemp * temprhs * increyxi 


u0 = np.r_[2.49282627e-12, 2.49282621e-12, 2.49282004e-12, 2.49220314e-12, 2.43131611e-12, 2.35213266e-13, 100.]
t = np.r_[1.0, 10.]
p = np.r_[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]


rhs_address = rhs.address

@njit
def time_func():
    usol, success = lsoda(rhs_address, u0, t, data=p, rtol=1.0e-3, atol=1.0e-6)
    assert success

time_func()
print(usol)