There is a small sign error in the trigonometric term in your solution for $\dot{r}$. A complete solution follows the sake future readers.
Problem statement
Is there a periodic solution for the following dynamical system?
$$
%
\begin{align}
%
\dot{x} &= 4 x+2 y - x\left(x^2+y^2\right)\\
%
\dot{y} &= -2 x+y-y \left(x^2+y^2\right)
%
\tag{1}
\end{align}
%
$$
Solution method
Use the theorem of Poincare and Bendixson to identify a trapping region, here the gray annulus where the sign of the radial time derivative can change.
The invariant region must
- Be closed and bounded,
- Not contain any critical points.
Solution
Identify critical points
At what points
$
\left[
\begin{array}{c}
x \\
y \\
\end{array}
\right]
$ does
$
\left[
\begin{array}{c}
\dot{x} \\
\dot{y} \\
\end{array}
\right]
=
\left[
\begin{array}{c}
0 \\
0 \\
\end{array}
\right]
$?
The only critical point is the origin.
Switch to polar coordinates
The workhorse formula is
$$
%
\begin{align}
%
x &= r \cos \theta, \\
%
y &= r \sin \theta.
%
\end{align}
%
$$
With $r^{2} = x^{2} + y^{2}$, use implicit differentiation to find
$$
r\dot{r} = x \dot{x} + y \dot{y}
\tag{2}
$$
Compute $\dot{r}$
Substituting into $(2)$ using $(1)$, and noting $\cos^{2} \theta = \frac{1}{2} \left( 1 + \cos 2\theta \right)$,
$$
%
\begin{align}
%
r \dot{r} &= x \dot{x} + y \dot{y} \\
%
&= r^{2} - r^{4} \color{blue}{+} 3x^{2} \\
%
&= r^{2} + \frac{3}{2} r^{2} \left( 1 \color{blue}{+} \cos 2\theta \right) - r^{4}
%
\end{align}
%
$$
Therefore
$$
\dot{r} = -r^{3} + \frac{r}{2} \left( 5 \color{blue}{+} 3 \cos 2\theta \right)
\tag{3}
$$
Classify $\dot{r}$
Look for regions where the flow is outward $\dot{r}>0$, and regions where the flow is inward $\dot{r}<0$.
(Note the interesting comment by @Evgeny.)
Classify the problem by examining the limiting cases of $\cos 2\theta$ at $\pm 1$
Outward flow: $\cos 2 \theta \ge -1$
$$
%
\begin{align}
%
\dot{r}_{in} &= -r^{3} + \frac{r}{2} \left( 5 + 3 (-1) \right) \\
%
&= r(1-r^{2})
%
\end{align}
%
$$
When $r<1$, $\dot{r}>0$, and the flow is outward.
Inward flow: $\cos 2 \theta \le 1$
$$
%
\begin{align}
%
\dot{r}_{in} &= -r^{3} + \frac{r}{2} \left( 5 + 3 (-1) \right) \\
%
&= r(4-r^{2})
%
\end{align}
%
$$
When $r>2$, $\dot{r}<0$, and the flow is inward.
Trapping region
The region between the two zones is the annulus centered at the origin with inner and outer radii
$$
%
\begin{align}
%
r_{in} &= 1 \\
%
r_{out} &= 2
%
\end{align}
%
\tag{4}
$$
There are no critical points. There region is closed and bounded. Therefore, a periodic solution exists.
Visualization
The vector field $\left[
\begin{array}{c}
\dot{x} \\
\dot{y} \\
\end{array}
\right]$ in $(1)$ is plotted against the gray trapping region in $(4)$. The red, dashed lines are nullclines which intercept at the critical point.
For the computation of the return or Poincaré map we can set $\theta=t$. Then the question to solve is if the integration of $r$ from $0$ to $2\pi$ can return to the same value $r_*$.
The radius equation is a Bernoulli equation, giving
$$
\frac{d}{dt}r^{-2}=-2r^{-3}\dot r=-2(r^{-2}-1+r^{-2}\cos t)
$$
This now is linear in $v=r^{-2}$. The integrating factor is $\phi(t)=e^{2t+2\mu\sin t}$,
$$
\frac{d}{dt}(\phi v)=2\phi, \\
\phi(2\pi)v(2\pi)-\phi(0)v(0)=2\int_0^{2\pi}e^{2t+2\mu\sin t}\,dt
$$
For a closed solution $v(0)=v(2\pi)=r_*^{-2}$, giving the equation for the fixed point
$$
r_*^{-2}=\frac2{e^{4\pi}-1}\int_0^{2\pi}e^{2t+2\mu\sin t}\,dt
$$
This will give large positive values $v_0=r_0^{-2}$ for large $\mu$, but has no singularities. Thus a limit cycle always exist, with increasingly small but always positive radius for increasing values of the $\mu$ parameter.
Best Answer
Consider a periodic orbit $p(t) = (x(t),y(t))$. By the Jordan curve theorem, $p$ splits the plane in two connected components, one of which (call it $D$) is bounded. Either $D$ is a subset of the annulus $A$ or $D$ contains the ball $B$ of radius 1 around the origin (think of a 'winding number' of the orbit being zero or nonzero that determines which case we are in). We will show that the first case cannot happen.
Define a vector field $F$ by $$ F(x,y) = \begin{pmatrix} -(x+y-y^3) \\ x-y-x^3 \end{pmatrix} = \begin{pmatrix} F_1(x,y) \\ F_2(x,y) \end{pmatrix}. $$ Note that $F(p(t)) \cdot \dot{p}(t) = -y'x'+x'y' = 0$. Stokes' theorem tells us $$ 0 = \int_0^T F(p(t)) \cdot \dot{p}(t) dt = \iint_D \left( \frac{\partial F_2}{\partial x} - \frac{\partial F_1}{\partial y} \right) dxdy = \iint_D ( 2 - 3 (x^2+y^2) ) dxdy. $$ The integrand $2 - 3 (x^2+y^2) =: g(x,y)$ is strictly negative in the annulus $A$, so this integral being zero implies that $D$ cannot be a subset of $A$.
Now suppose there were two distinct periodic orbits $p_1$ and $p_2$, which bound $D_1$ and $D_2$ respectively. Since both $D_1$ and $D_2$ contain the ball $B$, they are not disjoint. At the same time, the boundaries of $D_1$ and $D_2$ cannot intersect because distinct orbits cannot intersect. Thus, $D_1$ is a subset of $D_2$ or vice versa. But $$ \iint_{D_2 \backslash D_1} g(x,y) dxdy = \iint_{D_2} g(x,y) dxdy - \iint_{D_1} g(x,y) dxdy = 0-0 = 0 $$ and $D_2 \backslash D_1$ is a subset of $A$ on which $g$ is strictly negative, so $D_2 \backslash D_1$ has Lebesgue measure zero. But two periodic orbits cannot bound a set of measure zero between them (more precisely, $D_2 \backslash \bar{D}_1$ is nonempty and open, hence, has positive Lebesgue measure). This derives a contradiction to the existence of two distinct periodic orbits.
I do not recognise the theorem you quote about asymptotic stability of the orbit, but it adds up nicely as follows. I assume the $f$ in that theorem is the vector field generating the system, that is $f = (F_2,-F_1)$. Then $\mathrm{div}(f) = g$ is the integrand from before, which is strictly negative in the entire annulus $A$, so we can deduce that $\int_0^T \mathrm{div}(f)(p(t)) dt$ is strictly negative without explicitly knowing the orbit $p$.