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.
With the help of MATHEMATICA.
$$
A = \left(
\begin{array}{cccc}
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
-k & -k-K & 0 & 0 \\
k & -k & 0 & -\delta \\
\end{array}
\right)
$$
and the command
P = LyapunovSolve[A, Transpose[A], -IdentityMatrix[4]]
we obtain
$$
P = \left[
\begin{array}{cccc}
-\frac{\delta ^2 \left(3 k K+k (3 k+2)+K^2+K\right)+K (k+K) (2 k+K+1)}{2 \delta k (k+K) (2 k+K)} & -\frac{\delta ^2 k+\delta ^2 K+K}{4 \delta k^2+2
\delta k K} & -\frac{1}{2} & -\frac{k+1}{2 (k+K)} \\
-\frac{\delta ^2 k+\delta ^2 K+K}{4 \delta k^2+2 \delta k K} & \frac{\delta ^2 (k+K)+K (2 k+K+1)}{2 \delta (k+K) (2 k+K)} & \frac{k+1}{2 (k+K)} &
-\frac{1}{2} \\
-\frac{1}{2} & \frac{k+1}{2 (k+K)} & -\frac{\frac{\delta ^2 \left(2 k^2+2 k K+k+K^2\right)}{k+K}+(k+1) K}{2 \delta k} & \frac{K}{2 \delta } \\
-\frac{k+1}{2 (k+K)} & -\frac{1}{2} & \frac{K}{2 \delta } & \frac{k K+K}{2 \delta k+2 \delta K} \\
\end{array}
\right]
$$
Best Answer
There is no asymptotic stability. Taken the system
$$ \cases{ \dot x_1 = x_2\\ \dot x_2 = -x_1^3 } $$
and multiplying by $x_1^3, x_2$ as
$$ \cases{ x_1^3\dot x_1 = x_1^3x_2\\ x_2\dot x_2 = -x_1^3x_2 } $$
and adding we have
$$ \frac 12x_1^4+x_2^2 = C $$
so those orbits characterize a center around the origin.