Writes in polar coordinates $x''-(1-x^2-(x')^2)x'+4x=0$. I need to determine a periodic solution for :$x''-(1-x^2-(x')^2)x'+4x=0$.
We have the equivalent system:
$$\begin{cases} 
       x'=y \\
       y'=(1-x^2-y^2)y-4x.
   \end{cases}
$$

*

*We determined the stationary points for the equivalent system: $(0,0)$.


*I saw that this solution is unstable.


*But now I have to turn this system into polar coordinates. Unfortunately, I tried in all possibilities and I failed to bring the system to a beautiful shape depending on the polar coordinates.
I present what I tried:
Polar coordinates:
$\begin{cases} 
       x(r,\theta)=r\cos(\theta) \\
       y(r,\theta)=r\sin(\theta)
   \end{cases}
$
We have that $x^2+y^2=r^2$ and $\tan(\theta)=\frac{y}{x}$. So $r'=\frac{xx'+yy'}{r}$ and $\theta'=\frac{xy'-x'y}{r^2}.$
So $$r'=\frac{xy+y[(1-x^2-y^2)y-4x]}{r}=\frac{xy+y^2-x^2y^2-y^4-4xy}{r}.$$
I tried to calculate, but I have no idea how to bring the system into polar coordinates to continue my work.
Thanks!
 A: Note $$\theta'=\frac{(1-x^2-y^2)xy-4x^2-y^2}{r^2}=((1-r^2)\sin\theta-3\cos\theta)\cos\theta-1$$ and $$r'=\frac{xy+(1-x^2-y^2)y^2-4xy}r=-3r\sin\theta\cos\theta+(1-r^2)r\sin^2\theta$$ so $$(\log r)'=((1-r^2)\sin\theta-3\cos\theta)\sin\theta.$$ Thus $$\frac{(\log r)'}{\theta'+1}=\tan\theta.$$
A: The term in the middle acts as friction for $(x,x')$ outside the unit circle and as anti-friction inside the unit circle. This means that any particle moving under this dynamic will move away from the origin if close to it, and towards the unit circle if far away from it. The periodic solution will thus be close to the unit circle.
The numerical search for a periodic orbit can be formulated as boundary value problem, having the end point of the integration equal to the initial point. To account for the unknown period, map the whole problem to the constant interval $[0,1]$ using a parameter transformation.
In python/scipy this can be realized as follows
def odesys(t,u,p):
    x,v = u
    T, = p
    a = (1-x**2-v**2)*v-4*x
    return T*np.array([v,a])

def bc(u0, u1, p): return [u0[0], u1[0], u1[1]-u0[1]]

T = 2*np.pi
s = np.linspace(0,1,21);
u = [ np.sin(T*s), np.cos(T*s) ]

res = solve_bvp(odesys,bc,s,u,[T])
print(res.message, res.p)

s = np.linspace(0,1,301)
u = res.sol(s)
plt.plot(u[0],u[1],'b')
plt.plot(res.y[0],res.y[1],'xb')

This gives the result

