Circular trajectories

The system of ODEs $$\dot{x} = xy,\qquad \dot{y} = \frac{1-x^2+y^2}{2},$$ has circular trajectories, except those on the $y$-axis. Click your mouse on the plane, and the circular trajectories with the movement of the points will be shown.

The solution with intial condition $x(0)=x_0,y(0)=y_0$ is given by $$x(t) = \frac{2x_0}{1+x_0^2+y_0^2+(1-x_0^2-y_0^2)\cos t - 2y_0\sin t},\quad y(t) = \frac{(1-x_0^2-y_0^2)\sin t + 2y_0\cos t}{1+x_0^2+y_0^2+(1-x_0^2-y_0^2)\cos t - 2y_0\sin t}.$$ (for interested readers) One way to find this explicit solution is to add one new equation $\dot{z} = -yz$ with the initial condition $z(0)=1$ with the following steps:
1. Show that $xz$ is conserved, to deduce that $xz=x_0$.
2. Substitute $x=x_0/z$ and $y = -\dot{z}/z$ into the equation $\dot{y}=\frac{1-x^2+y^2}{2}$ to get a second order ODE for $z$.
3. From the second order ODE for $z$, show that $\frac{\dot{z}^2+z^2+x_0^2}{z}$ is invaraint, and hence $\frac{\dot{z}^2+z^2+x_0^2}{z} = 1+x_0^2+y_0^2.$ Here you need the equation $\dot{z} = -yz$ to get $\dot{z}(0)=-y_0$.
4. The solution of the above first order ODE $\frac{\dot{z}^2+z^2+x_0^2}{z} = 1+x_0^2+y_0^2$ is $z(t) = \frac{1+x_0^2+y_0^2}{2} + A\cos t + B\sin t$, for some constants $A$ and $B$. Use the initial condition for $z(0)$ and $\dot{z}(0)$ to find the solution $z(t)$, and then $x(t),y(t)$.