There is only one maximal Lyapunov exponent (MLE). In my answer I will tell you how to calculate this, and not many Lyapunov exponents (because it is a different method).
The method I am presenting is due to Benettin [1], and simply evolves in parallel 2 trajectories while constantly rescaling one of the two.
MLE calculation
Two neighboring trajectories with initial distance d0 are evolved in time.
At time $t_i$ their distance $d(t_i)$ either exceeds either some upper_threshold,
or is lower than some lower_threshold, which initializes
a rescaling of the test trajectory back to having distance d0 from
the given one, while the rescaling keeps the difference vector along the maximal
expansion/contraction direction:
$$
u_2 = u_1 + (u_2 - u_1)/(d(t_i)/d_0)
$$
where $u_2, u_1$ are the test and reference states respectively (both vectors).
The maximum
Lyapunov exponent is the average of the time-local Lyapunov exponents
$$
\lambda = \frac{1}{t_{n}}\sum_{i=1}^{n}
\ln\left( a_i \right),\quad a_i = \frac{d(t_{i})}{d_0}.
$$
Software that does it
DynamicalSystems.jl is a open-sourced software library that offers functions that can calculate:
- maximum lyapunov exponent.
- lyapunov spectrum (all Lyapunov exponents). The library also describes how this method works in detail, in case your question was about computing many Lyapunov exponents instead of only the maximum one. The method you describe about how to find the MLE of a 1D map can be expanded into the method described in the link.
- MLE from a numerical dataset.
[1] : G. Benettin et al., Phys. Rev. A 14, pp 2338 (1976)