clear
clf
c=input(' Give c value as complex number: ');
kmax=input(' Give maximum applications of inverse mapping: ');
c2=sqrt(1-4*c);
if(real(c2)>=0)
  zunst=(1+c2)/2;
else
  zunst=(1-c2)/2;
end
z=zunst;
for i=1:kmax
  root=sqrt(z-c);
  rn=rand;
  if(rn<0.5)
    z=root;
  else
    z=-root;
  end
  x(i)=real(z);
  y(i)=imag(z);
end
plot(x,y,'.','markersize',3);
axis([-2 2 -2 2])
axis square