Lyapunov Functions and Basins of Attraction \begin{align}
    x' &= -x^3 + x^5 + (x^4)(y^5)\\[.7em]
    y' &= -8y^3 + y^5 - 10(y^4)(x^5)
\end{align}
$(0,0)$ is obviously a critical point of the system, and we are given that it is asymptotically stable, but have to show it. 
I have tried to make a Lyapunov function $V(x,y) = ax^2 + cy^2$, with a,c > 0 but I am having trouble to prove that $\frac{d}{dt} V(x,y)$ is negative definite. I get some complicated polynomial I can't use logic to finalize. How can I change the Lyapunov to come up with a meaningful conclusion?
\begin{align}
    \frac{d}{dt}V(x,y) = 2ax(-3x^2 + 5x^4 + 4x^3y^5) + 2cy(-24y^2+5y^4-40y^3x^5)\\[.7em]
\end{align}
 A: Hint: choose $a$ and $c$ so the terms in $x^5 y^5$ cancel.  Note that for $x$ near $0$, lower powers of $x$ dominate higher powers, and similarly for $y$.
A: For 
$$V(x,y)=\frac{1}{2}(x^2+y^2)$$
we have
$$\dot{V}=-x^4-8y^4+x^6+y^6-9y^5x^5\\
\leq -(1-x^2-|xy|)x^4-8(1-(1/8)y^2-|xy|)y^4$$
Thus if we define the region (neighborhood of the origin)
$$\Omega:=\left\{(x,y)|(x^2+|xy|<1)\textrm{ and }\frac{1}{8}y^2+|xy|<1\right\}$$
we have that $$\dot{V}(x,y)<0\qquad \forall (x,y)\in\Omega\neq 0$$
Now if we choose the level set 
$$\Omega_0:=\left\{(x,y)|V(x,y)<\frac{1}{4}\right\}\subseteq \Omega$$
we also have that 
$$\dot{V}(x,y)<0\qquad \forall (x,y)\in\Omega_0\neq 0$$
Therefore every solution that starts within $\Omega_0$ remains therein and asymptotically converges to the origin ($V$ is strictly decreasing).
A: This problem can be handled with an optimization procedure, having in mind that generally is a non convex problem. The result depends on the test Lyapunov function used so we will generalize to a quadratic Lyapunov function
$$
V(p) = p^{\dagger}\cdot M\cdot p = a x^2+b x y + c y^2,\ \ \ p = (x,y)^{\dagger}
$$
and
$$
f(p) = \{-x^3 + x^5 + x^4 y^5, -8 y^3 + y^5 - 10 x^5 y^4\}
$$
with $a>0,c>0, a b-b^2 > 0$ to assure positivity on $M$. We will assure a set involving the origin $Q_{\dot V}$ such that $\dot V(Q_{\dot V}) < 0$. The optimization process will be used to guarantee a maximal $Q_{\dot V}$.
After determination of $\dot V = 2 p^{\dagger}\cdot M\cdot f(p)$ we follow with a change of variables
$$
\cases{
x = r\cos\theta\\
y = r\sin\theta
}
$$
so $\dot V = \dot L(a,b,c,r,\theta)$. The next step is to make a sweep on $\theta$ calculating
$$
S(a,b,c, r)=\{\dot V(a,b,c,r,k\Delta\theta\},\ \ k = 0,\cdots, \frac{2\pi}{\Delta\theta}
$$
and then the optimization formulation follows as
$$
\max_{a,b,c,r}r\ \ \ \ \text{s. t.}\ \ \ \ a > 0, c> 0, a c -b^2 > 0, \max S(a,b,c,r) \le -\gamma
$$
with $\gamma > 0$ a margin control number.
Follows a MATHEMATICA script which implements this procedure in the present case.
f = {-x^3 + x^5 + x^4 y^5, -8 y^3 + y^5 - 10 x^5 y^4};
V = a x^2 + 2 b x y + c y^2;
dV = Grad[V, {x, y}].f /. {x -> r Cos[t], y -> r Sin[t]};
rest = Max[Table[dV, {t, -Pi, Pi, Pi/30}]] < -0.2;
rests = Join[{rest}, {r > 0, a > 0, c > 0, a c - b^2 > 0}];
sols = NMinimize[Join[{-r}, rests], {a, b, c, r}, Method -> "DifferentialEvolution"]
rest /. sols[[2]]

dV0 = Grad[V, p].f /. sols[[2]]
V0 = V /. sols[[2]]
r0 = 1.5;
rmax = r /. sols[[2]];
gr0 = StreamPlot[f, {x, -r0, r0}, {y, -r0, r0}];
gr1a = ContourPlot[dV0, {x, -r0, r0}, {y, -r0, r0}, ContourShading -> None, Contours -> 80];
gr1b = ContourPlot[dV0 == 0, {x, -r0, r0}, {y, -r0, r0}, ContourStyle -> Blue];
gr2 = ContourPlot[x^2 + y^2 == rmax^2, {x, -r0, r0}, {y, -r0, r0}, ContourStyle -> {Red, Dashed}];
Show[gr0, gr1a, gr1b, gr2]

Follows a plot showing in black the level sets $Q_{\dot V}$ an in blue the trace of $\dot V = 0$. In dashed red is shown the largest circular set $\delta = 0.99992$ defining the maximum attraction basin for the given test Lyapunov function's family.

