clear
kmax=input(' Give maximum applications of mapping: ');
x=-2.2:.01:2; % Create grid of initial c values
y=x;
[xx,yy]=meshgrid(x,y);
s=length(x);
k=0;
x2=xx.^2;
y2=yy.^2;
c=xx+i*yy;
while(k<=kmax)
   yy=2*xx.*yy+imag(c); % Do mapping
   xx=x2-y2+real(c);
   k=k+1;
   x2=xx.^2;
   y2=yy.^2;
end
[igrid,jgrid]=find( x2+y2 <= 4); % Find indices of points still in R after kmax
xplot=x(jgrid)';
yplot=y(igrid)';
plot(xplot,yplot,'r.','markersize',4)
axis([-2 2 -2 2])
axis square