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.
Best Answer
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. $$