Prove that the system has a single limit cycle We have the following system:
\begin{cases} 
      x'=x(x^2+y^2-2x-3)-y\\
      y'=y(x^2+y^2-2x-3)+x
   \end{cases}
We want to determine that it has a single limit cycle.
To prove this we had to use the Poincare Bendixson theorem.
We have the linearized system:
\begin{cases}
   x'=-y\\
   y'=x
\end{cases}
$A=\begin{pmatrix}
0 & -1 \\
1 & 0 
\end{pmatrix}$. We can see that we have a center point.
After that I found the stationary points. My problem arose when I had to determine the set A, which does not contain stationary points and is invariant to flow. I tried to solve this problem, but without any success.
 A: Caution: In its present state, this "hint" provides only a part of the answer. I am going to revise it in a short while.
Hint: consider dot product
$$(x,y).(x',y')=xx'+yy'=(x^2+y^2)(x^2+y^2-2x-3)$$
It is negative inside the circle with equation $x^2+y^2-2x-3=0$, i.e. with center $(1;0)$ and radius $2$,
This geometric property, means that a particle submitted to the corresponding field is "directed" in the half plane (defined by the tangent) where the circle's center is situated. Only a partial conclusion can be derived here.

Matlab program:
function edh
clear all;close all;hold on;axis equal off;
axis([-1.5,4.5,-3,3]);
t=0:pi/100:2*pi;
plot(1+2*exp(i*t),'c','linewidth',2);% circle
plot([-1.5,4.5],[0,0]);plot([0,0],[-3,3]);
tspan = 0:0.01:10;  % temporal span
ode = @(t,y) sc(y); 
[x,y] = meshgrid(-1.5:0.25:4.5,-3:0.25:3);
w=x.^2+y.^2-2*x-3;% a way to describe the diff. system
u=x.*w-y;
v=y.*w+x;
quiver(x,y,u,v,2);% arrows
for r=1.8:0.1:2.2
    for k=-pi:pi/20:pi
        y0 = [1+r*cos(k);r*sin(k)]; % initial points
        [t,y] = ode45(ode, tspan, y0); % Runge Kutta    
        plot(y(:,1),y(:,2),'color','r');
    end
end
%%%
function yp = sc(y)
w=y(1).^2+y(2).^2-2*y(1)-3; % another way to 
u=y(1).*w-y(2); % describe the diff. system
v=y(2).*w+y(1);
yp=[u;v];

A: As was established, for the dynamic of the radius you get
$$
\dot r=r(r^2-2x-3).
$$
This can be weakened to a scalar differential inequality
$$
r((r-1)^2-4)\le\dot r\le r((r+1)^2-4)\\
r(r-3)(r+1)\le r\le r(r-1)(r+3).
$$
The lower bound is positive for $r>3$, the upper bound negative for $r<1$. This means that the annulus $1\le r\le 3$ is a trapping region in negative time direction. As the angular velocity around the origin is always $1$, there are no equilibrium points except the origin. This establishes the existence of a repelling limit cycle inside the annulus. However not that it is single, there could in theory also be a triple of repelling, attracting and again repelling limit cycles.

Numerical experiments suggest that the limit cycle is close to the circle of radius $2$ around $(1,0)$, that is the root locus of the factor $x^2+y^2-2x-3$. Let $R=\sqrt{(x-1)^2+y^2}$, then
$$\begin{align}
\dot RR&=(x-1)\dot x+y\dot y=((x-1)x+y^2)(R^2-4)-(x-1)y+yx\\
&=R^2(R^2-4)+(x-1)(R^2-4)+y
\end{align}$$
By Cauchy-Schwarz
$$
|(x-1)(R^2-4)+y|\le R\sqrt{(R^2-4)^2+1}
$$
so that
$$
R(R^2-4)-\sqrt{(R^2-4)^2+1}\le \dot R\le \underbrace{R(R^2-4)}_{=f(R)}+\underbrace{\sqrt{(R^2-4)^2+1}}_{=g(R)}
$$

The upper bound is negative on $[1.1,1.8]$, the lower bound positive for $R>2.2$, so that the annulus $1.8\le R\le 2.2$ is a trapping region in negative time direction. The plot of the numerical solution shows that this bound is rather tight.

