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.
You could pick the control law as follows
$$
u(x) = -x_1 - x_2 - \frac{x_1^2}{x_2},
$$
but that is not well defined when $x_2=0$.
Instead one could make use of the fact that your proposed control law makes the system dynamics linear, with
$$
\dot{x} = A\,x,
$$
$$
A = \begin{bmatrix}0 & 1 \\ -1 & -1\end{bmatrix}.
$$
A Lyapunov function for such system can be found of the form
$$
V(x) = x^\top P\,x,
$$
with $P$ positive definite which satisfies the Lyapunov equation
$$
A^\top P + P\,A = -Q,
$$
with $Q$ positive definite. If $A$ is stable any positive definite $Q$ should also yield a corresponding positive definite $P$.
For example when setting $Q$ equal to the identity matrix yields
$$
V(x) = \frac{1}{2} (3\,x_1^2+2\,x_1\,x_2+2\,x_2^2) = \frac{1}{2} (2\,x_1^2+(x_1+x_2)^2+x_2^2),
$$
$$
\dot{V}(x) = -x_1^2 - x_2^2.
$$
Best Answer
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.