Formal proof of Lyapunov stability I was trying to solve the question of AeT. on the (local) Lyapunov stability of the origin (non-hyperbolic equilibrium) for the dynamical system
$$\dot{x}=-4y+x^2,\\\dot{y}=4x+y^2.\tag{1}$$
The streamplot below indicates that this actually is true.

Performing the change of variables to polar coordinates $x=r\cos\phi$, $y=r\sin\phi$ and after some trigonometric manipulations we result in
$$\dot{r}=r^2(\cos^3\phi+\sin^3\phi)\\ \dot{\phi}=4+r\cos \phi \sin\phi(\cos \phi -\sin \phi )$$
From this set of equations I want to prove that if we start with sufficiently small $r$ then $r$ will remain bounded with very small variations over time.
My intuitive approach: For very small $r$
$$\dot{\phi}\approx 4$$ that yields $$\phi(t)\approx 4t +\phi_0$$
If we replace in the $r$ dynamics we obtain
$$\dot{r}\approx r^2\left[\cos^3(4t+\phi_0)+\sin^3(4t+\phi_0)\right]$$
Integrating over $[0,t]$ we obtain
$$\frac{1}{r_0}-\frac{1}{r(t)}\approx \int_0^t{\left[\cos^3(4s+\phi_0)+\sin^3(4s+\phi_0)\right]ds}$$
The right hand side is a bounded function of time with absolute value bounded by $4\pi$ since
$$\int_{t_0}^{t_0+2\pi}{\left[\cos^3(4s+\phi_0)+\sin^3(4s+\phi_0)\right]ds}=0 \quad \forall t_0$$
Thus for very small $r_0$ it holds true that $r(t)\approx r_0$.
I understand that the above analysis is at least incomplete (if not erroneous) and I would be glad if someone can provide a rigorous treatment on the problem.
I think that a "singular-perturbation like" approach may be the solution (bounding $r$ by $\epsilon$) and considering the comparison system to prove the global boundedness result but I haven't progressed much up to now.

*

*OP's streamplot suggests that the line $y=x-4$ is a flow trajectory.  If we insert the line $y=x-4$ in OP's eq. (1) we easily confirm that this is indeed the case.

*From now on we will assume that $y\neq x-4$. It is straightforward to check that the function
$$H(x,y)~:=~\frac{xy+16}{x-y-4}-4 \ln |x-y-4| $$
is a first integral/an integral of motion: $\dot{H}=0$.

*In fact, if we introduce the (non-canonical) Poisson bracket
$$B~:=~\{x,y\}~:=~ (x-y-4)^2 ,$$
then OP's eq. (1) becomes Hamilton's equations
$$ \dot{x}~=~\{x,H\}, \qquad \dot{y}~=~\{y,H\}. $$

*The above result was found by following the playbook laid out in my Phys.SE answer here: $B$ is an integrating factor for the existence of the Hamiltonian $H$.
