Consider the system for $(x_1(t),x_2(t))$
\begin{align} \dot{x}_1 &= x_1^2+x_1^3+x_2\\ \dot{x}_2 &= x_1^2+u \end{align}
Find a function $u=\phi(x_1(t),x_2(t))$ so that if
$$V(x) = \frac12(x_1^2+(x_1+x_1^2+x_1^3+x_2)^2)$$
then $\dot{V}=-2V$.
So my method was to find $\dot V$ then equate that to $-2V$ and rearrange to find $u$. But this seems very messy, is there a nicer way to do this?
\begin{align} \dot{V} &= x_1\dot{x}_1 + (\dot x_1 + 2x_1\dot x_1 + 3x_1^2\dot x_1 + \dot x_2)(x_1+x_1^2+x_1^3+x_2) \\ &= \dot x_1\big[x_1 + (1 + 2x_1+3x_1^2)(x_1+x_1^2+x_1^3+x_2)\big] + \dot x_2(x_1+x_1^2+x_1^3+x_2) \\ &= (x_1^2+x_1^3+x_2)\big[x_1 + (1 + 2x_1+3x_1^2)(x_1+x_1^2+x_1^3+x_2)\big] + (x_1^2+u)(x_1+x_1^2+x_1^3+x_2) \\ &= (x_1^2+x_1^3+x_2)\big[x_1 + (1 + 2x_1+3x_1^2)(x_1+x_1^2+x_1^3+x_2)\big] + x_1^2(x_1+x_1^2+x_1^3+x_2) + u(x_1+x_1^2+x_1^3+x_2)\\ &=-(x_1^2+(x_1+x_1^2+x_1^3+x_2)^2) \end{align}
$$\iff u(x_1+x_1^2+x_1^3+x_2)=-(x_1^2+(x_1+x_1^2+x_1^3+x_2)^2)-x_1^2(x_1+x_1^2+x_1^3+x_2)-(x_1^2+x_1^3+x_2)\big[x_1 + (1 + 2x_1+3x_1^2)(x_1+x_1^2+x_1^3+x_2)\big]$$
$$\implies u=-\frac{(x_1^2+(x_1+x_1^2+x_1^3+x_2)^2)}{(x_1+x_1^2+x_1^3+x_2)}-x_1^2-(x_1^2+x_1^3+x_2)\left[\frac{x_1}{(x_1+x_1^2+x_1^3+x_2)} + (1 + 2x_1+3x_1^2)\right]$$
This is clearly an example of backstepping design. This is a sequential design methodology for strict feedback nonlinear systems (note that in your example $x_2$ can be seen as an 'input' to the $x_1$ subsystem). So the idea is to select an appropriate virtual control $\alpha_1(x_1)$ that can stabilize the $x_1$ subsystem if $x_2=\alpha_1(x_1)$ and then drive the difference $z_2:=x_2-\alpha_1(x_1)$ to zero by selecting suitable $u$.
So if you select $$\alpha_1(x_1)=-x_1^2-x_1^3-x_1$$ and define the new error variable $$z_2:=x_2-\alpha_1(x_1)=x_1+x_1^2+x_1^3+x_2$$ then your $V$ takes the form $$V=\frac{1}{2}x_1^2+\frac{1}{2}z_2^2$$ and $$\dot{x}_1=-x_1+z_2\\ \dot{z}_2=x_1^2+u-\frac{d\alpha_1}{dx_1}\dot{x}_1$$ Note now that $$\dot{V}=x_1(-x_1+z_2)+z_2\dot{z}_2=-x_1^2+z_2(x_1+x_1^2+u-\frac{d\alpha_1}{dx_1}\dot{x}_1)$$ Thus selecting $$u=-x_1-x_1^2+\frac{d\alpha_1}{dx_1}\dot{x}_1-z_2$$ guarantees $$\dot{V}=-x_1^2-z_2^2=-2V$$ Hence the desired control input is $$u=-x_1-x_1^2+(-1-2x_1-3x_1^2)(x_1^2+x_1^3+x_2)-(x_2+x_1+x_1^2+x_1^3)$$