I’ve come across a problem regarding using lyapunovs() with stiff diffeq solvers. It is very similar to this one:

Problem with tangent_integrator: Chaostools tangent_integrator stiff solvers

https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/87

Since lyapunovs() also uses a tangent_integrator, the same problem arises when one wants to calculate the lyapunov exponents with a stiff solver. These solvers seem to be incompatible with the tangent_integrator. Turning autodiff off switches to numerical differentiation, which is super slow. Can anyone help me with a workaround (the link to the workaround in the issue above is no longer available), maybe leaving the variables untyped, so it doesn’t run into problems when converting Float64 into Dual, idk ,or tips to make numerical differentiation faster (Jacobi matrix optimalization etc)?

Minimal Working Example:

```
function lorenz(du,u,p,t)
sigma = p[1]
rho = p[2]
beta = p[3]
# eom
du[1] = sigma*(u[2]-u[1])
du[2] = u[1]*(rho-u[3])-u[2]
du[3] = u[1]*u[2]-beta*u[3]
end
function lorenz_jac(J,u,p,t)
sigma = p[1]
rho = p[2]
beta = p[3]
J[1,1] = -sigma
J[1,2] = sigma
J[1,3] = 0
J[2,1] = rho-u[3]
J[2,2] = -1
J[2,3] = -u[1]
J[3,1] = u[2]
J[3,2] = u[1]
J[3,3] = -beta
end
state = rand(3)
param = [16.,45.92,4.]
ds = ContinuousDynamicalSystem(lorenz,state,param,lorenz_jac)
lyapunovs(ds,2000.;Ttr = 2000.,alg = Rosenbrock23())
```

In this case, it is possible to calculate the lyapunov exponents in a matter of minutes using numerical differentiation (autodiff = false). However in my case, the Jacobian looks something like this:

```
function njacobi(J,u,p)
a = -p[1]*p[2]/(p[4]-1)
c = -p[1]*(1-p[2])
n = convert(Int32,p[4])
x = u[1:2:end]
y = u[2:2:end]
for i in 1:(2*n)
for j in 1:(2*n)
if i<(n+1) && j<(n+1)
if i!=j
J[i,j] = a*x[i]
else
J[i,j] = p[1]*(1 - (1 - p[2])*y[i] - p[2]/(n-1)*skipsum_jacobi(u,i))
end
elseif i<(n+1) && j>n
if (j-i) == n
J[i,j] = c*x[i]
else
J[i,j] = 0
end
elseif i>n && j<(n+1)
if (i-j) == n
J[i,j] = p[3]
else
J[i,j] = 0
end
else
if i==j
J[i,j] = -p[3]
else
J[i,j] = 0
end
end
end
end
end
```

which is supposed to look like this:

```
4×4 Array{Float64,2}:
4.92501 -2.27912 -5.7178 0.0
-1.36462 4.99815 0.0 -3.42351
0.1 0.0 -0.1 0.0
0.0 0.1 0.0 -0.1
```

for a 4D version of my system. It can be treated as 4 2X2 block matrices, and the elements are given depending on which block I’m in, hence the conditionals in my Jacobian. The two bottom blocks are always constant and diagonal. The loops are there because I want to change the dimensions from time to time using a dimension parameter in the parameter container. The problem is, in this case, calculating lyapunovs via numerical differentiation runs for ages, especially when I try to increase the dimensions (50,100+). How can this be done properly?