Let's start out using notation similar to the example you linked to:
$$ \ddot{y}+b\dot{y}+\sin(y)=a\cos(ct)+d $$
As in the example, we'll write this in autonomous form:
$$ \mathbf{x}=\begin{pmatrix}x_1\x_2\x_3\end{pmatrix}=\begin{pmatrix}y\\dot{y}\ct\end{pmatrix} $$
the differential equation can now be written:
$$ \mathbf{\dot{x}}=\begin{pmatrix}\dot{x_1}\\dot{x_2}\\dot{x_3}\end{pmatrix}=\mathbf{f(x)}=\begin{pmatrix}x_2\-bx_2-\sin{x_1}+a\cos{x_3}+d\c\end{pmatrix} $$
Now, suppose we study the evolution of two different starting points in this system that are arbitraily close. Let's denote the trajectory of one of them with $\mathbf{x}(t)$ and the other one with $\mathbf{x}(t)+\mathbf{z}(t)$, where $\mathbf{z}(t)$ is the vector difference between the two points. The time derivative of $\mathbf{z}(t)$ can be written:
$$ \mathbf{\dot{z}}=\mathbf{f(x+z)}-\mathbf{f(x)} $$
If $\mathbf{z}$ is sufficiently small this can be written:
$$ \mathbf{\dot{z}}=\mathbf{Jz} $$
where $\mathbf{J(x)}$ is the Jacobian of $\mathbf{f(x)}$:
$$ \mathbf{J}=\begin{pmatrix} 0 & 1 & 0\ -\cos{x_1} & -b & -a\sin{x_3}\ 0 & 0 & 0 \end{pmatrix} $$
Over a short period $\delta t$, the change in $\mathbf{z}$ is given by:
$$ \mathbf{z}(t+\delta t)=\mathbf{z}(t)+\delta t\mathbf{Jz}(t)=(\mathbf{I} +\delta t\mathbf{J})\mathbf{z}(t) $$
Assuming $\mathbf{z}$ and $\delta t$ are small enough for the linearizations to hold we can decribe the change in $\mathbf{z}$ from time $t_1$ to $t_2$ with:
$$ \mathbf{z}(t_2)=\left[\prod_i(\mathbf{I} +\delta t_i\mathbf{J}(\mathbf{x}_i))\right]\mathbf{z}(t_1) $$
where $\delta t_i$ are sufficiently small time intervals from $t_1$ to $t_2$ and $\mathbf{x}_i$ is $\mathbf{x}$ at those times.
What you need to do to calculate the Lyaponov exponents of the this system is to numerically solve the trajectory $\mathbf{x}(t)$ for an arbitrary starting point $\mathbf{x}_0$ over a long time interval $T$. Then you need to calculate the matrix $\mathbf{A}$:
$$ \mathbf{A}=\prod_i(\mathbf{I} +\delta t_i\mathbf{J}(\mathbf{x}_i)) $$
using the values for $\delta t_i$ and $\mathbf{x}_i$ from the solved trajectory. 804
You then need to calculate the eigenvalues of $\mathbf{A}$. Let's call them $\alpha_1$, $\alpha_2$ and $\alpha_3$. The lyaponov exponents are given by:
$$ \lambda_i=\frac{1}{T}\log{\alpha_i} $$
assuming $T$ is long enough. You need to experiment with different $\mathbf{x}_0$ and $T$. Hopefully, the variations in the Lyaponov exponents you calculate for different $\mathbf{x}_0$ will decrease with larger $T$ and converge for sufficiently large $T$. 915