Lyapunov Exponent of the Rössler System

Question

I would like to plot a semi-log plot of the Lyapunov exponent $λ_{x_0}, x_0=(1,1,1) \in \mathbb{R}^3$ versus time $t$, for the Rössler system described by: \begin{equation} \begin{array}{lcl}\dot{x}& =& -y-z \\ \dot{y} &=&x+ay \\ \dot{z}&=& b+z(x-c) \end{array} \end{equation} where $a,b,c \in \mathbb{R}$ real parameters of the system.

Now, the attractor is usually met for $a=b=0.2$ and $c=5.7$. My target is to show that the plot which is derived by plotting the log of the distance between trajectories is well approximated by a straight line with a positive slope (which is actually the Lyapunov exponent), so the Rössler system has a positive Lyapunov exponent and experiences chaotic motion.

What I have managed so far by searching in stack.exchange is to find something similar on Lyapunov Exponent for Rössler System which is provided as:

ClearAll["Global`*"];
deq1 = -(y1[t] + z1[t]);
deq2 = x1[t] + 0.1 y1[t];
deq3 = 0.2 + x1[t] z1[t] - 5.7 z1[t];

deq4 = -(y2[t] + z2[t]);
deq5 = x2[t] + 0.1 y2[t];
deq6 = 0.2 + x2[t] z2[t] - 5.7 z2[t];

x10 = 1; y10 = 1; z10 = 1;
dx0 = 10^-8;
x20 = x10 + dx0; y20 = y10; z20 = z10;
tin = 0; tfin = 10000;
tstep = 1;
acc = 12;

lcedata = {};
sum = 0;

d0 = Sqrt[(x10 - x20)^2 + (y10 - y20)^2 + (z10 - z20)^2 ];

For[i = 1, i < tfin/tstep, i++,

sdeq = {x1'[t] == deq1, y1'[t] == deq2, z1'[t] == deq3, 
x2'[t] == deq4, y2'[t] == deq5, z2'[t] == deq6, x1[0] == x10, 
y1[0] == y10, z1[0] == z10, x2[0] == x20, y2[0] == y20, 
z2[0] == z20};
sol = NDSolve[
sdeq, {x1[t], y1[t], z1[t], x2[t], y2[t], z2[t]}, {t, 0, tstep}, 
MaxSteps -> Infinity, Method -> "Adams", PrecisionGoal -> acc, 
AccuracyGoal -> acc];

xx1[t_]  = x1[t] /. sol[[1]];
yy1[t_]  = y1[t] /. sol[[1]];
zz1[t_] = z1[t] /. sol[[1]];

xx2[t_]  = x2[t] /. sol[[1]];
yy2[t_]  = y2[t] /. sol[[1]];
zz2[t_] = z2[t] /. sol[[1]];

d1 = Sqrt[(xx1[tstep] - xx2[tstep])^2 + (yy1[tstep] - yy2[tstep])^2 + 
          (zz1[tstep] - zz2[tstep])^2 ];

sum += Log[d1/d0];
dlce = sum/(tstep*i);
AppendTo[lcedata, {tstep*i, Log10[dlce]}];

w1 = (xx1[tstep] - xx2[tstep])*(d0/d1); 
w2 = (yy1[tstep] - yy2[tstep])*(d0/d1);
w3 = (zz1[tstep] - zz2[tstep])*(d0/d1); 

x10 = xx1[tstep];
y10 = yy1[tstep];
z10 = zz1[tstep];

x20 = x10 + w1;
y20 = y10 + w2;
z20 = z10 + w3;

i = i++;

If[Mod[tstep*i, 100] == 0, 
Print[" For t = ", tstep*i, " , ", " LCE = ", dlce]]

]

S0 = ListPlot[{lcedata}, Frame -> True, Axes -> False, 
PlotRange -> All, Joined -> True, 
          FrameLabel -> {"t", "log10(LCE)"}, 
FrameStyle -> Directive["Helvetica", 17], ImageSize -> 550]

But what I would like to show is something similar to this analysis Lyapunov exponent of the Lorenz system.Am I doing something wrong?

Thanks!