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
So you found $$ \frac{d}{dt}(x_1^2+x_2^2)=-2x_2^2-2x_2^4=-x_2^2(1+2x_2^2)-\dot x_1x_2 $$ Now apply integration-by-parts to the last term to get $$ \frac{d}{dt}(x_1^2+x_2^2+x_1x_2) =-x_2^2(1+2x_2^2)+x_1\dot x_2 \\=-x_2^2(1+2x_2^2)-x_1^2-x_1x_2(1+x_2^2) \\\le-x_1^2-x_2^2(1+2x_2^2)+\frac{x_1^2+x_2^2}2(1+x_2^2) $$ which is negative outside the origin for $|x_2|<1$. Note that the modified quadratic form is still positive definite.
Remark: If you take the derivative of the second equation and insert the first one you get $$ \ddot x_2+(1+3x_2^2)\dot x_2+x_2=0 $$ which is a Liénard equation and the trick above is standard there.