Getting unexpected values for Characteristic Lyapunov Exponents in Python - What could be going wrong with my implementation of Benettin algorithm?

I wrote a code to calculate the Characteristic Lyapunov Exponents in Python for a Lorenz system. I for sure did something wrong because the CLEs are way off but I cannot figure out what. I am using the Benettin et al. algorithm. I start by randomly choosing and initial state, and then 3 random tangent vectors that I orthonormalize using Gram Schmidt process (which also calculates the volumes). After the initialisation, I start iterating. I evolve the reference trajectory using Euler method and the tangent vectors using w(t+1)=L@w(t) where L is the Jacobian matrix at x(t) and w(t) is the tangent vector. I orthonormalize every northo and calculate the sum of log(volumes). Every iwrite the CLEs are printed out. If anyone has any idea what could be the issue, I'd be grateful. I am implementing this algorithm for a Lorenz system to validate it in order to use it for a different system.

Here is the code:

import numpy as np

def lorenz(x, y, z, sigma, rho, beta):
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

def lorenz_L(x,y,z):
    return np.array([
    [-sigma, sigma, 0],
    [rho - z, -1, -x],
    [y, x, -beta]])   

def gram_schmidt(V):
    """
    Orthogonalize a set of vectors stored as the ROWS of matrix V.
    It also computes the volume generated the first i vectors for i=1,...,dim
    """
    # Get the number of vectors.
    n = V.shape[0]
    W = np.copy(V)
    R = np.zeros(n)
    VOL = np.zeros(n)
    for j in range(n):
        # To orthogonalize the vector in row j with respect to the
        # previous vectors, subtract from it its projection onto
        # each of the previous vectors.
        for k in range(j):
            W[j] -= np.dot(W[k], W[j]) * W[k]
        R[j] = np.linalg.norm(W[j])
        W[j] = W[j] / R[j]
        VOL[0]=R[0]
    for k in range(1,n):
        VOL[k]=VOL[k-1]*R[k]
    return W, VOL

# Parameters
sigma = 16.0
rho = 45.92
beta = 4.0

t_start = 0.0
t_end = 100.0
steps = 10000
northo = 1
iwrite = 1000 
num_points = steps * northo
dt = t_end/num_points
icont=0
dim=3

X = np.random.random(dim)

# Estimate Lyapunov exponents

num_lyapunov_exponents = dim
lyapunov_exponents = np.zeros(num_lyapunov_exponents)

num_vec=dim

sum_volumes = np.zeros(dim)
sum_lyapunov_exponents = np.zeros(dim)
cle = np.zeros(dim)

tangent_vec = np.random.random((dim,dim))
tangent_vec, _ = gram_schmidt(tangent_vec) 

for j in range(steps):
    for is_ in range(northo):
        icont=icont+1
        # velocities
        dX_dt = lorenz(X[0],X[1],X[2],sigma,rho,beta)

        # evolution of the reference trajectory
        for i in range(dim):
            X[i] = X[i] + dX_dt[i] * dt
 
        # stability matrix for Lorenz system
        L = lorenz_L(X[0],X[1],X[2])
        
        # evolve tangent vectors
        tangent_vec = (L@tangent_vec.T).T
        # for ivec in range(num_vec):
        #     tangent_vec[ivec] = L@tangent_vec[ivec]
    
    #orthonormalize
    tangent_vec, volumes = gram_schmidt(tangent_vec) 

    for i in range(dim):
        sum_volumes[i]=sum_volumes[i]+np.log(volumes[i])
    
    # calculate and print results
    if icont%iwrite == 0:
        t = icont*dt
        print(f'time : {t}, dt : {dt}, frame : {icont}')
        sum_lyapunov_exponents = sum_volumes / t
        cle[0]=sum_lyapunov_exponents[0]
        for l in range(1,dim):
            cle[l]=sum_lyapunov_exponents[l]-sum_lyapunov_exponents[l-1]
        for ii in range(dim):
            print(f'{ii+1}   {cle[ii]}')
        print()

I tried the above code, and I get Lyapunov exponents depending on what I set the time stepping values. I expect to get the theoretical expectation of a Kaplan York dimension of 2.07, First lyapunov exponent positive (1.5;2), the second zero and the third negative (-20;-30).