function dy=phase_focus_fun(t,y)

%Stable
mu=1;
%Unstable
%mu=-1;

dy1=-y(2)-mu*y(1)*(y(1)^2+y(2)^2);
dy2=y(1)-mu*y(2)*(y(1)^2+y(2)^2);

dy=[dy1; dy2];