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.
The theorem works: as @Ian wrote in his comment, the flux across the lower surface must be negative. The right computation is $$\int \int_{S_3}(x+y,y,z^2) \cdot n dS = \int_0^{2 \pi} \int_0^1(*,*,1) \cdot(0,0,-\rho)d\rho d\theta = - \pi$$ since $S_3$ has radius $1$. Notice that the outward normal is pointing out, while yours is probably pointing inside the volume $V$ defined by your surface.
The other flux across $S_1$ is $16 \pi$, which can be done as the one above. That said, your net flux using the definition is $\frac{14}{3} \pi - \frac{15}{2} \pi + 16 \pi - \pi = \frac{14}{3} \pi + \frac{15}{2}\pi$
There's another problem in your divergence: that integral is $$\int \int \int_V 2 + 2z dxdydz = \int_1^2 \int_{B_z(0,0)}2+2z dxdydz =[\ldots ] = \frac{14}{3} \pi +\frac{15}{2} \pi$$
which confirms the previous result.
Best Answer
It is better to note that $$ \frac12(x^2+y^2)'=xx'+yy'=x^2+y^2-(x^4+y^4) $$ is positive for very small $\|(x,y)\|$ and negative for very large $\|(x,y)\|$.