To my knowledge there is not a general method for finding a Lyapunov function. In this case one can solve the differential equations and use that to find a Lyapunov function. Namely $x_2$ is decoupled from $x_1$ and can be shown to have the following solution
$$
x_2(t) = C_1\,e^{-t},
$$
where $C_1$ is a constant and depends on the initial condition of $x_2$. Substituting the above equation into the expression for $\dot{x}_1$ gives
$$
\dot{x}_1 = x_1 (C_1\,e^{-t} -1)
$$
which is a separable differential equation, namely
$$
\frac{dx_1}{x_1} = (C_1\,e^{-t} -1) dt.
$$
Integrating on both sides gives
$$
\log(x_1) = -C_1\,e^{-t} -t+C_2.
$$
Solving for $x_1$ gives
\begin{align}
x_1(t) &= e^{-C_1\,e^{-t} -t+C_2}, \\ &= C_3\,e^{-C_1\,e^{-t} -t}, \\
&= C_3\,e^{-t}\,e^{-C_1\,e^{-t}},
\end{align}
or when using the definition for $x_2$ then it can also be expressed as $x_1(t)=C_3\,e^{-t}\,e^{-x_2}$. So the quantities $x_2$ and $x_1\,e^{x_2}$ will both decay exponentially fast, so the following Lyapunov function can be used
$$
V(x) = x_2^2 + x_1^2\,e^{2\,x_2},
$$
for which it can be shown that its derivative is
$$
\dot{V}(x) = -2\,x_2^2 - 2\,x_1^2\,e^{2\,x_2}.
$$
I will leave proving that $V(x)$ is radially unbounded to you.
These are multiple question, I will try to answer some of them.
I assume that $a, b > 0$.
An equilibrium is isolated if there does not exist a small "ball" around it such that it is the only equilibrium in that ball. For example:
$$
\dot{x_1} = 0
$$
has no isolated equilibrium because no matter how close you look at an equilibrium, there are always infinitely many other equilibrium points close by. Similiar, the system
$$
\begin{align}
\dot{x}_1 &= x_1 - x_2 \\
\dot{x}_2 &= x_1 - x_2 \\
\end{align}
$$
has no isolated equilibrium: For example $x_1 = x_2 = 1$ is is an equilibrium, but so is $x_1 = x_2 = 1.1$ and $x_1 = x_2 = 1.01$ and $x_1 = x_2 = 1.00000001$ and so on. No matter how close you look, there are always "more" equilibrium points.
So to your example, you have:
$$
\begin{align}
\dot{x}_1 &= x_2 + x_3 \\
\dot{x}_2 &= -a \sin(x_1) - b x_2 \\
\dot{x}_3 &= -a \sin(x_1) + x_3
\end{align}
$$
To get the equilibrium points you set $\dot{x} = 0$ so you need $x_2 = -x_3$ from the first equation. Then you get
$$
\begin{align}
\dot{x}_2 &= \dot{x}_3 = 0 \\
\implies -a \sin(x_1) - b x_2 &= -a \sin(x_1) + x_3 \\
\implies -b x_2 &= x_3 \\
\implies -b x_2 &= -x_2
\end{align}
$$
We assumed that $b > 0$ and if $b \neq 1$ the last equation can only be true if $x_2 = 0$ and then $x_3 = 0$. If you insert that you get
$$
-a \sin(x_1) = 0
$$
so the equilbrium points are $x_1 = \sin(k \pi), k \in \mathbb{Z}, x_2 = 0, x_3 = 0$. Although that means you have infinitly many equilibrium points they are all isolated because the roots of $\sin(x_1) = 0$ are isolated. Conclusion: You can use Lyapunov method.
If $b = 1$ the system has no isolated equilibrium points.
You can use indirect method but it will only tell you local results and sometimes you can't use it (if you get zero eigenvalues).
You have the function
$$
V(x) = \frac{1}{2} x^T x + a(1 - \cos(x_1))
$$
but it is not a Lyapunov function for your system:
$$
\dot{V}(x) = - b x_2^2 + x_1 x_2 - x_3^2 + x_1 x_3
$$
If $x_3 = 0$ and $x_1 = 2 b x_2$ you get $\dot{V}(x) = b x_2^2$ and that is positive no matter how small you make $x_2$. So this function can't be used.
Best Answer
From the Jacobian that you obtained for the stable equilibrium points $(-1, 0)$, and $(1, 0)$, the linearized system becomes
$$ \left[\matrix{\dot{x}_{1} \cr \dot{x}_{2}}\right] = \left[\matrix{0 & 1 \cr -2 & -\delta}\right] \left[\matrix{x_{1} \cr x_{2}}\right] .$$
Consider the Lyapunov function in the following quadratic form:
$$ V(\mathbf{x}) = \frac{1}{2} \mathbf{x}^{T} \mathbf{P} \mathbf{x} = \frac{1}{2} \left[\matrix{x_{1} & x_{2}}\right] \left[\matrix{\frac{\delta^{2} + 2 \delta + 2}{\delta} & 1 \cr 1 & \frac{\delta + 1}{\delta}}\right] \left[\matrix{x_{1} \cr x_{2}}\right] .$$
The positive definite $\mathbf{P}$ is found by solving the Lyapunov equation $\mathbf{P} \mathbf{A} + \mathbf{A}^{T} \mathbf{P} = - \mathbf{Q}$. In this case, $\mathbf{Q}$ is chosen as
$$ \mathbf{Q} = \left[\matrix{2 & 0 \cr 0 & \delta}\right] .$$
The derivative of $V(\mathbf{x})$ along the trajectories of the linear system is given by
$$ \dot{V}(\mathbf{x}) = - 2 x_{1}^{2} - \delta x_{2}^{2} < 0 \quad \text{for} \quad \mathbf{x} \neq \mathbf{0} .$$
Remember that this result, however, is local.
Region of Attraction
To estimate the region of attraction, the previous Lyapunov function can be used to determine a domain $D$ about the stable equilibrium by solving the derivative of $V(\mathbf{x})$ along the trajectories of the nonlinear system (Duffing Oscillator) where $\dot{V}(\mathbf{x})$ is negative definite.
$$ \dot{V} = \left(\frac{\delta^{2} + 2 \delta + 2}{\delta}\right) x_{1} \dot{x}_{1} + \dot{x}_{1} x_{2} + x_{1} \dot{x}_{2} + \left(\frac{\delta + 1}{\delta}\right) x_{2} \dot{x}_{2} $$
$$ \dot{V} = \left(\frac{\delta^{2} + 2 \delta + 2}{\delta}\right) x_{1} x_{2} + x_{2}^{2} + x_{1} \left(x_{1} - x_{1}^{3} - \delta x_{2}\right) + \left(\frac{\delta + 1}{\delta}\right) x_{2} \left(x_{1} - x_{1}^{3} - \delta x_{2}\right) $$
$$ \dot{V} = - x_{1}^{4} - \left[\left(\frac{\delta + 1}{\delta}\right) x_{1} x_{2} - 1\right] x_{1}^{2} + \left(\frac{3 \delta + 3}{\delta}\right) x_{1} x_{2} - \delta x_{2}^{2} $$
For example, if $\delta = 1$, the region of attraction is found to be satisfied in $D_{(1, 0)} = \left\{\sqrt{2} < x_{1} < \sqrt{5}, x_{2} \in \mathbb{R}\right\}$, and $D_{(-1, 0)} = \left\{-\sqrt{5} < x_{1} < -\sqrt{2}, x_{2} \in \mathbb{R}\right\}$.
The streamplot for the Duffing Oscillator with $\delta = 1$ is shown below. Note that $\pm \sqrt{5} \approx \pm 2.236$.