n = 60;
dfi = 2*pi/n;
Fi = dfi*(0:n)';
c = ['yrgb'];
for nr = 1:4
   a = 0.8 * nr;
   fi = 0;
   R = 1; r= 1;
   for j = 1:n
       fi = fi + dfi;
       h = 1;
       while abs(h) > 1e-10
       f = r^2 - 1-a*sin(0.5*r*r*sin(2*fi));
       fprim = 2*r-a*r*cos(0.5*r*r*sin(2*fi))*sin(2*fi);
       h = -f/fprim;
       r = r+h;
       end
       R = [R;r];
   end
      x = R.*cos(Fi); y = R.*sin(Fi); plot(x,y,c(nr)), hold on
    
end
axis equal











