exclude errors when computing the Lyapunov exponents.

In short, we here remind that the trajectories can be classified as function of the Lyapunov exponents. Thus, one can distinguish periodical behaviours, for $\lambda=0$, dissipative movements with a fixed point or a basin of attraction $\left(\lambda<0\right)$, and aperiodical chaotic states $\left(\lambda>0\right)$, when the iterative discrete evolution of the solution series (Eqs. 5 and 6) leads to a chaotic pattern.

Another way of measuring the system sensitivity to initial conditions is to compute the largest Lyapunov exponent ([LLE]{}). Usually a couple of methods can be employed, one based on the time dependence of the multi-dimensional Lyapunov exponent, the other on Wolf’s standard method that uses a Gram-Schmidt Reorthonormalization of the tangent vectors [@wolf-85]. In the latter, the [LLE]{} is obtained by taking the asymptotic value of the multi-dimensional Lyapunov exponent:

$$\lambda_{Max}=\lim_{t \rightarrow \infty} \frac{\lambda\left(t\right)}{t}.$$

Still, this method has the disadvantage that the integration times have to be at least an order of magnitude larger than those here considered. Other methods are slightly less efficacious, being more CPU time-consuming when simulating strong chaotic systems [@ramasubramanian-00].

We consequently used the first method and noticed the saturation behaviour, *i.e.* the arising of a plateau after a 