xmin=-2;
xmax=1;
ymin=-1.5;
ymax=1.5;

npts=1000;
cx=xmin:(xmax-xmin)/npts:xmax;
cy=ymin:(ymax-ymin)/npts:ymax;
[CX,CY]=meshgrid(cx,cy);

C = CX + i*CY;
Z=C;
cnt=ones(size(C));
jmax=100;

init=cputime;
for j=1:jmax
    Z = Z.^2 + C;
    inside=abs(Z)<=2;
    cnt=cnt+inside;
end
cnt=log(cnt);
fin=cputime-init
pcolor(CX,CY,cnt),shading interp
colormap(flipud(jet(jmax)))
%colormap(jet)
axis square
axis off
