N = 500; M = 500;
y = [];

z1 = 1;
z2 = 1/2*(-1+sqrt(-3));
z3 = 1/2*(-1-sqrt(-3));

for i = 1:N, i
for j = 1:M,

zstart = 2*(i/N-1/2) + sqrt(-1)*2*(j/M-1/2);
z = zstart;
for idx = 1:100,

	znew = 1/(3*z^2) + 2/3 * z;	
	if abs(znew - z) < 1/200,z = znew;break; end;
	z = znew;
end;	

if abs(z-z1) < 1/100,
	y(i,j) = 1;
elseif abs(z-z2) < 1/100,
	y(i,j) = 2;
elseif abs(z-z3) < 1/100,
	 y(i,j) = 3;	
else
	disp('Error'); zstart
end;	

end;
end;

imagesc(y',[1,3]); axis('off');