# Problem with calculation of Lyapunov exponent

The code given by Chris K for Lyapunov Exponent does not work and gives lot of errors for the dynamical system given in equation(2) of this paper

The dynamical equations:

eqns = {x'[t] == 5*(y[t] - x[t]) + 15*y[t]*z[t],
y'[t] == -10*(y[t]*Abs[y[t]]) - y[t] + 4*x[t]*z[t],
z'[t] == 3.9*z[t] - x[t] y[t]};


with ics:

ics = {x -> 0, y -> 5, z -> 0};


The phase space diagram of the attractor:

The borrowed code for the Lyapunov exponent is given below:

GramSchmidt[w_?MatrixQ] :=
Module[{v = ConstantArray[0, Length[w]]},
Table[v[[n]] =
w[[n]] -
Sum[(v[[i]] . w[[n]]/v[[i]] . v[[i]])*v[[i]], {i, n - 1}], {n,
Length[w]}];
v]

LyapunovExponents[eqnsin_List, icsin : ({__Rule} | _Association),
nlein_Integer : 0, opts___?OptionQ] :=
Module[{(*options*)tstep, maxsteps, ndsolveopts, logbase, showplot,
plotexponents, plotopts,(*iterators*)c, i,
j,(*other variables*)\[Delta], neq, nle, vars, rhs, jac, eqns,
unks, ics, cum, res, edat, state, newstate, sol, W,
norms},(*parse options*)
tstep = Evaluate[
TStep /. Flatten[{opts, Options[LyapunovExponents]}]];
maxsteps =
Evaluate[
MaxSteps /. Flatten[{opts, Options[LyapunovExponents]}]];
ndsolveopts =
Evaluate[
NDSolveOpts /. Flatten[{opts, Options[LyapunovExponents]}]];
logbase =
Evaluate[LogBase /. Flatten[{opts, Options[LyapunovExponents]}]];
showplot =
Evaluate[
ShowPlot /. Flatten[{opts, Options[LyapunovExponents]}]];
plotexponents =
Evaluate[
PlotExponents /. Flatten[{opts, Options[LyapunovExponents]}]];
plotopts =
Evaluate[
PlotOpts /. Flatten[{opts, Options[LyapunovExponents]}]];
neq = Length[eqnsin];
If[nlein == 0, nle = neq,
nle = nlein];(*how many exponents*)(*extract vars and right hand \
sides from eqnsin*)vars = eqnsin[[All, 1, 0, 1]];
rhs = eqnsin[[All, 2]];
(*jacobian matrix*)jac = D[rhs, {Replace[vars, {x_ -> x[t]}, 1]}];
eqns =
Join[eqnsin,
Flatten[Table[\[Delta][i, j]'[
t] == (jac . Table[\[Delta][i, j][t], {i, neq}])[[i]], {j,
nle}, {i, neq}]]];
unks =
Join[vars, Flatten[Table[\[Delta][i, j], {j, nle}, {i, neq}]]];
ics = Join[Table[var[0] == (var /. icsin), {var, vars}],
Flatten[Table[\[Delta][i, j][0] ==
IdentityMatrix[neq][[i, j]], {j, nle}, {i, neq}]]];
cum = Table[0, {nle}];
state =
First@NDSolveProcessEquations[Flatten[Join[eqns, ics]], unks, t,
Evaluate[Sequence @@ ndsolveopts]];
(*main loop*)
edat = Table[newstate = First@NDSolveReinitialize[state, ics];
NDSolveIterate[newstate, c tstep];
sol = NDSolveProcessSolutions[newstate];
W = GramSchmidt[
Evaluate[
Table[\[Delta][i, j][c tstep], {j, nle}, {i, neq}] /. sol]];
norms = Map[Norm, W];
(*update running vector magnitudes*)
cum = cum + Log[logbase, norms];
ics =
Join[Table[var[c tstep] == (var[c tstep] /. sol), {var, vars}],
Flatten[Table[\[Delta][i, j][c tstep] == (W/norms)[[j, i]], {j,
nle}, {i, neq}]]];
cum/(c tstep), {c, maxsteps}];
If[showplot,
Print[ListPlot[Transpose[edat][[1 ;; plotexponents]],
Evaluate[Sequence @@ plotopts]]]];
Return[cum/(maxsteps tstep)]];

Options[LyapunovExponents] = {NDSolveOpts -> {}, TStep -> 1,
MaxSteps -> 10^4, LogBase -> E, ShowPlot -> False,
PlotExponents -> 3, PlotOpts -> {}};

LyapunovExponents[eqns, ics, ShowPlot -> True]


Errors:

How to resolve this type of error and get the desire result?

In your differential equations, replace Abs with RealAbs.