Hello there,

I’m trying to figure out how to solve my lyapunov spectrum handwrite code using the QR method for a simple example (lorenz system), there is something wrong with my code ? I’m trying to reproduce the values obtained in this topic : https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent

Here is my code:

```
using DynamicalSystems
using Plots, DataFrames, DifferentialEquations, LinearAlgebra
# Lorenz system
function loop(u, p, t)
σ = p[1]; ρ = p[2]; β = p[3]
du1 = σ*(u[2]-u[1])
du2 = u[1]*(ρ-u[3]) - u[2]
du3 = u[1]*u[2] - β*u[3]
return SVector{3}(du1, du2, du3)
end
# Jacobian:
function loop_jac(u, p)
σ, ρ, β = p
J = @SMatrix [-σ σ 0;
ρ - u[3] (-1) (-u[1]);
u[2] u[1] -β]
return J
end
#initial parameters
u0 = [1,1,1]
p0 = [10,28,8/3]
t_span = (0,100.0)
prob = ODEProblem(loop, u0, t_span,p0)
solution = solve(prob, RK4(), reltol=1e-11, abstol=1e-11, dense=true)
global U = [1 0 0; 0 1 0;0 0 1]
df = DataFrame(solution)
lyapuA=[]
lyapuB=[]
lyapuC=[]
append!(lyapuA,1)
append!(lyapuB,1)
append!(lyapuC,1)
AA = 0
BB = 0
CC = 0
for i in 2:length(df.timestamp)
v0 =solution.u[i]
dt = df.timestamp[i] - df.timestamp[i-1]
U_n = (U + loop_jac(v0,p0)*dt)*U
Q, R = qr(U_n)
U = Q
AA = log(abs(diag(R)[1]))/(df.timestamp[i])
BB = log(abs(diag(R)[2]))/(df.timestamp[i])
CC = log(abs(diag(R)[3]))/(df.timestamp[i])
append!(lyapuA,AA)
append!(lyapuB,BB)
append!(lyapuC,CC)
end
p2 = plot(df.timestamp,lyapuA)
p3 = plot(df.timestamp,lyapuB)
p4 = plot(df.timestamp,lyapuC)
plot(p2,p3,p4)
```