c=0.4+0.07*i;
x=0;
for k=1:7
   x=x^2+c;
end
rho=abs(x);
clear x
t=linspace(0,2*pi,500);
X=rho*(cos(t)+i*sin(t));
plot(real(X),imag(X),'k.');
hold on
%Two preimages
for k1=1:2
   Z1=((-1)^k1)*sqrt(X-c);
   plot(real(Z1),imag(Z1),'b.');
   for k2=1:2
      Z2=((-1)^k2)*sqrt(Z1-c);
      plot(real(Z2),imag(Z2),'g.');
   end
end
plot(real(c),imag(c),'r*');
hold off
clear Z1 Z2
pause
plot(real(X),imag(X),'k.');
hold on

%Four preimages
for k1=1:2
   Z1=((-1)^k1)*sqrt(X-c);
   plot(real(Z1),imag(Z1),'b.');
   for k2=1:2
      Z2=((-1)^k2)*sqrt(Z1-c);
      plot(real(Z2),imag(Z2),'g.');
      for k3=1:2
         Z3=((-1)^k3)*sqrt(Z2-c);
         plot(real(Z3),imag(Z3),'y.');
         for k4=1:2
            Z4=((-1)^k4)*sqrt(Z3-c);
            plot(real(Z4),imag(Z4),'m.');
        end
      end
   end
end
plot(real(c),imag(c),'r*');
hold off
clear Z1 Z2 Z3 Z4
pause

plot(real(X),imag(X),'k.');
hold on
for k1=1:2
   Z1=((-1)^k1)*sqrt(X-c);
   plot(real(Z1),imag(Z1),'b.');
   for k2=1:2
      Z2=((-1)^k2)*sqrt(Z1-c);
      plot(real(Z2),imag(Z2),'g.');
      for k3=1:2
         Z3=((-1)^k3)*sqrt(Z2-c);
         plot(real(Z3),imag(Z3),'y.');
         for k4=1:2
            Z4=((-1)^k4)*sqrt(Z3-c);
            plot(real(Z4),imag(Z4),'m.');
            for k5=1:2
               Z5=((-1)^k5)*sqrt(Z4-c);
               plot(real(Z5),imag(Z5),'k.');
               for k6=1:2
                  Z6=((-1)^k6)*sqrt(Z5-c);
                  plot(real(Z6),imag(Z6),'b.');
                  for k7=1:2
                     Z7=((-1)^k7)*sqrt(Z6-c);
                     plot(real(Z7),imag(Z7),'g.');
                  end
               end
            end
         end
      end
   end
end
plot(real(c),imag(c),'r*');
hold off