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.
Well, if the trajectories are confined to some region which does contain a fixed point, then it may happen that they all converge to that fixed point, so that there are no limit cycles in the region. So by assuming that there are no fixed points, you rule out that possibility.
But note that when you apply Poincaré–Bendixson, the trapping region is usually an annulus (or something topologically equivalent to that), so you don't consider the whole region inside the potential limit cycle. As you say, there is necessarily a fixed point inside the inner circle of the annulus, but that point doesn't belong to the annulus.
Best Answer
Using the change of variables $(u,v)=(x+y,y)$, the $(x,y)$-differential system is equivalent to the $(u,v)$-differential system $$u'=uv^2-2v,\qquad v'=v^3+u-2v$$ In particular, $$(u^2+2v^2)'=2(uu'+2vv')=2v^2(u^2+2v^2-4)$$ This shows that the ellipsis $(E)$ of equation $u^2+2v^2=4$ is invariant by the dynamics. In the $(x,y)$-plane, the equation of $(E)$ is $$x^2+2xy+3y^2=4.$$ Starting from every point inside $(E)$, one converges to $(x_\infty,y_\infty)=(0,0)$. Starting from every point outside $(E)$, one diverges in the sense that $x(t)^2+y(t)^2\to+\infty$. Finally, starting from every point on $(E)$, one cycles on $(E)$ counterclockwise with time period $$\int_0^{2\pi}\frac{\mathrm dt}{\sqrt2-\sin(2t)}=2\pi.$$ To evaluate the period, recall that, for every $a\gt1$, $$\int_0^{2\pi}\frac{\mathrm dt}{a+\sin(t)}=\frac{2\pi}{\sqrt{a^2-1}}.$$
Edit: Some simulations of the systems $x'=x(y^2−c)−y$, $y'=y(y^2−c)+x$ for some $c\gt0$. The cycle of each such system is the ellipsis $(E_c)$ with equation $$x^2+2cxy+(1+2c^2)y^2=2c(1+c^2).$$ For $c=4$:
streamplot[{x(y^2-4)-y,y(y^2-4)+x},{x,-20,+20},{y,-5,+5}]
$\qquad\qquad$
For $c=.2$:
streamplot[{x(y^2-.2)-y,y(y^2-.2)+x},{x,-2,+2},{y,-2,+2}]
$\qquad\qquad$