# Lyapunov Exponent as a function of parameter

I am working at Rossler system and I would like to ask you for some help. This system is described by equations: \begin{array}{ll} \dot{x} = -y - z \\ \dot{y} = x + 0.2y \\ \label{eq:UkladRosSP} \dot{z} = 0.2 + z(x-c) \end{array} I am trying to plot a graph, on which a maximum Lyapunov Exponent would be shown. I have followed the way which was suggested in file : http://www.mathematica-journal.com/issue/v6i3/article/sandri/contents/63sandri.pdf . The problem is that the results does not converge to what it should be. For example, for c=1.1, we should get a 1-period cycle, while the MLE is non-negative (that means the chaos). Also the procedure of plotting takes a lot of time. Is there any way to get a better results? Code:

 RKStep[f_, y_, y0_, dt_] :=
Module[{k1, k2, k3, k4},
k1 = dt N[f /. Thread[y -> y0] ];
k2 = dt N[f /. Thread[y -> y0 + k1/2] ];
k3 = dt N[f /. Thread[y -> y0 + k2/2] ];
k4 = dt N[f /. Thread[y -> y0 + k3] ];
y0 + (k1 + 2*k2 + 2*k3 + k4)/6 ]

IntVarEq[F_List, DPhi_List, x_List, Phi_List, x0_List,
Phi0_List, {t1_, dt_}] :=
Module[{n, f, y, y0, yt},
n = Length[x0];
f = Flatten[Join[F, DPhi]];
y = Flatten[Join[x, Phi]];
y0 = Flatten[Join[x0, Phi0]];
yt = Nest[RKStep[f, y, #, N[dt]] &,
N[y0], Round[N[t1/dt]] ];
{First[#], Rest[#]} & @Partition[yt, n] ]

JacobianMatrix[funs_List, vars_List] := Outer[D, funs, vars]

Normi[x_] := Sqrt[x.x]

Lyap[c_] := Module[{},
F[{x_, y_, z_}] := {-z - y, x + 0.2*y, 0.2 + z*(x - c)};
x0 = {-1, 0, 0};
T = 40;
stepsize = 0.01;
n = Length[x0];
x = Array[a, n];
Phi = Array[b, {n, n}];
DPhi = Phi.Transpose[JacobianMatrix[F[x], x]];
Phi0 = IdentityMatrix[n];
{xT, PhiT} = IntVarEq[F[x], DPhi, x, Phi, x0, Phi0, {T, stepsize}];
xT
]


My aim is to achieve a plot similar to this one :

Any ideas?