% routine to create a Newton's method plot
% a VERY slow program, indeed

n = 21;   % an odd number prevents a zero
n = input('enter resolution ->');
dx = 3/n;
tol = 0.1;		% how close we need to be a root to
				% declare convergence


% w1, w2, w3 are the roots of the equation x^3-1=0

w1 = 1;
w2 = exp(2*pi*i/3);
w3 = w2*w2;

range = [-1.5 : dx : 1.5];


% mm records which root we're nearest

mm = zeros(n+1,n+1); 

for p=1:n+1, for q=1:n+1
	z = range(p) + i*range(q);

	while(1)	% do repeatedly:
		z = z - (z*z*z-1)/(3*z*z); % basic newton step
		if abs(z-w1) < tol
			mm(p,q)=1;
			break;
		end;
		if abs(z-w2) < tol
			mm(p,q)=2;
			break;
		end;
		if abs(z-w3) < tol
			mm(p,q)=3;
			break;
		end;
	end;
end;end;

% produce the image on the screen

image(mm');
axis('image');
axis('off');

% uncomment exactly one of these lines:

colormap(gray(3)); 	% for grayscale monitors
%colormap(eye(3));	% for color monitors