% routine to create a Newton's method plot

n  = 81;   % an odd number prevents a zero
n = input('enter resolution ->');
dx = 3/n;
nits = 10;		% number of full matrix iterations
tol = 0.25;		% 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;


disp('Generating base matrix');
[X,Y] = meshgrid(-1.5:dx:1.5);
m = X + i*Y;
disp(size(m));

clear X,Y;

disp(['Base matrix of size ',int2str(prod(size(m))),' created']);
disp(['Doing ',int2str(nits),' full matrix steps']);


% first round of calculations
% we compute the newton iteration for
% all entries of the matrix

for k=1:nits
	m = m - (m.^3-1)./(3*m.^2);
	disp(['Step #',int2str(k),' finished']);
end;

% mm records which root we're nearest

mm = 1*(abs(m-w1)<tol) + 2*(abs(m-w2)<tol) + 3*(abs(m-w3)<tol);

% the errs hold which entries are not converged yet

errs = find(mm==0);

nerr = prod(size(errs));
nels = prod(size(m));
disp(['Number of errors is ',int2str(nerr)]);
disp(['Error rate is ',num2str(nerr*100/nels),' %']);

xtra = 0;

% now we finish computing those entries which have yet to converge

[jj,kk] = size(mm);
for j = errs'
	x = m(j);
	xx = 1*(abs(x-w1)<tol)+2*(abs(x-w2)<tol)+3*(abs(x-w3)<tol);
	while xx==0
		x = x - (x*x*x-1) / (3*x*x);
		xtra = xtra+1;
		xx=1*(abs(x-w1)<tol)+2*(abs(x-w2)<tol)+3*(abs(x-w3)<tol);
	end; % end while
	m(j) = x;
	mm(j) = xx;
end;end;end;    % end for loops and if mm==0
disp([int2str(xtra),' extra Newton steps required']);

% produce the image on the screen

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

% uncomment exactly one of these lines:

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