Limit cycle of dynamical system $x'=xy^2-x-y$, $y'=y^3+x-y$ Consider a planar ODE system $z'=f(z)$ with $z=(x,y)$, 
$$
f(x,y)=(xy^2-x-y,y^3+x-y).
$$
Using polar coordinates, one can get 
$$
r'=r(r^2\sin^2\theta-1),\quad \theta'=1.
$$
With Mathematica one can get

As one can see from the figure, there is a limit cycle for the system. I tried to apply the Poincaré–Bendixson theorem to show the existence of the limit cycle by using some negatively invariant annular region. But there seems no hope of finding such region. How should I go on?

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$
