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,1000);
X=rho*(cos(t)+i*sin(t));
Z1=sqrt(X);
plot(real(X),imag(X),'b.');
hold on
plot(real(Z1),imag(Z1),'g.',real(-Z1),imag(-Z1),'g.');
Z11=sqrt(Z1);
Z12=sqrt(-Z1);
plot(real(Z11),imag(Z11),'r.',real(-Z11),imag(-Z11),'r.');
plot(real(Z12),imag(Z12),'r.',real(-Z12),imag(-Z12),'r.');