Stability for Nonlinear System I am trying to assess the (Liapunov) stability of the equilibrium at $(0,0)$ of the system
\begin{align*}
x_1' &= -4x_2 + x_1^2 \\
x_2' &= 4x_1 + x_2^2.
\end{align*}
I plotted the phase portrait in Mathematica, and it looks like a stable (but not asymptotically stable) equilibrium, with orbits circling about the origin.
However, my professor tells me the origin is actually unstable.
I don't know how to prove or disprove this rigorously. The equilibrium is non-hyperbolic, and I haven't been able to find a Liapunov function for the system. Any guidance would be appreciated. 
Edit: I tried polar coordinates: I get 
$$ r' = \frac1r(x_1x_1' + x_2x_2') = r^2(\cos^3 \theta + \sin^3 \theta)$$
and 
$$ \theta' = \frac{x_1x_2'-x_1'x_2}{r^2} = 4 + r\cos\theta\sin\theta(\sin\theta-\cos\theta).$$
So there's a half-plane where $r$ is increasing and a half-plane where it's decreasing. The angle $\theta$ is increasing in a neighborhood of the origin.
Is it reasonable to say that since 
$$ 0 = \int_{\theta_0}^{\theta_0 + 2\pi} \cos^3 \theta + \sin^3 \theta d\theta = \int_{\theta_0}^{\theta_0 + 2\pi} \frac{r'}{r^2} d\theta = \int_{\theta_0}^{\theta_0 + 2\pi} \left(\frac{1}{r}\right)' d\theta, $$
that the orbits are in fact periodic?
 A: Dynamical system
$$
%
\begin{align}
%
 \dot{x} &= -4y + x^{2} \\
%
 \dot{y} &= 4x + y^{2}
%
\end{align}
%
$$
The sole critical point is the origin.
Phase portrait
The nullclines are the dashed lines. Red for $\dot{x} = 0$, purple for $\dot{y} = 0$.

Analysis
Polar coordinate transformation produces
$$
  \dot{r} = r^{2} 
\left(
  \cos^{3} \theta + \sin^{3} \theta
\right)
$$ 
The plot following shows
$$
  f(\theta) = \cos^{3} \theta + \sin^{3} \theta
$$

The red zone shows contraction where $\color{red}{\dot{r} < 0}$, and the blue zone dilation, $\color{blue}{\dot{r} > 0}$

Notice that as we approach the origin, the dynamical system 
$$
%
\begin{align}
%
 \dot{x} &= -4y + x^{2} \qquad \Rightarrow \qquad -4y\\
%
 \dot{y} &= 4x + y^{2} \qquad \Rightarrow \qquad 4x
%
\end{align}
%
$$
approaches a periodic system.
A: This is not really an answer to the question. But I think that it can be useful for answering it.
If we divide the differential equations we obtain:
$$\dfrac{dx_2}{dx_1}=\dfrac{4x_1+x_2^2}{-4x_2+x_1^2}$$
In order to simplify notation I will replace $x_1 \to t$ and $x_2 \to x$
$$\dfrac{dx}{dt}=\dfrac{4t+x^2}{-4x+t^2}$$.
This differential equation does have a symmetry, with an infinitesimal generator
$$X = \dfrac{t-x-4}{t^2-4x}\partial_{x}.$$
By using this symmetry we can use the method of canonical coordinates to obtain the substitution:
$$t(r,s) = r$$
$$x(r,s)=r-4+\exp\left[W\left(-0.25(r^2-4r+16)\exp(0.25s(r)) \right)-0.25s(r) \right],$$
in which $W$ is the Lambert-$W$-Function.
Using this set of substitutions for the differential equation will result in 
$$\dfrac{ds(r)}{dr}=1\implies s(r) = r+c_1=t+c_1$$
and 
$$x(t)=t-4+\exp\left[W\left(-0.25(t^2-4t+16)\exp(0.25(t+c_1) \right)-0.25(t-c_1) \right].$$
Turning this back into the original coordinates we obtain:
$$x_2=x_1-4+\exp\left[W\left(-0.25(x_1^2-4x_1+16)\exp(0.25(x_1+c_1) \right)-0.25(x_1-c_1) \right]$$
A: Solving
$$
\frac{dx_2}{dx_1}=\frac{4x_1+x_2^2}{-4x_2-x_1^2}
$$
according to MATHEMATICA we get
$$
36 c_1 \left(x_1-x_2\right)-4 \left(36 c_1+2^{2/3} x_2\right)+2^{2/3} x_1 \left(x_2+4\right)+
   2^{2/3+2} \left(4-x_1+x_2\right) \ln \left(\left(x_1-4\right) x_1+16\right)+2^{2/3+3}
   \left(x_1-x_2-4\right) \tanh ^{-1}\left(\frac{\left(x_1-8\right) x_1+4 \left(x_2+8\right)}{x_1^2-4
   x_2}\right)=0
$$
and we can obtain a level set partition as shown in the following plot

typical of a conservative system.
