function  w = newtbas3(poly,wind,newtmap,count)

% Uses  count  iterations to draw Newton basins  
% in the viewing window  wind  for the polynomial  
% with descending vector  poly  of coefficients, 
% using the RGB color matrix  colors.  Try 
%   
%      poly = [1 -3 0 1];
%      wind = [-2 4 -2.25 2.25];
%      newtmap = [0.0  0.0  0.0;   	% black
%                 0.9  0.9  0.3;	% yellow
%                 0.8  0.0  0.1; 	% red
%                 0.2  0.7  0.3];	% green
%      count = 50;
%
% for the cubic polynomial  x^3 - 3x^2 + 1.

% Henry Edwards June 1998

r = roots(poly);	 %  Find polynomial roots
s = size(r);          n = s(1);
r1 = r(1);  r2 = r(2);  r3 = r(3);

a = wind(1);       b = wind(2);
c = wind(3);       d = wind(4);

a = wind(1);       b = wind(2);  % viewing window
c = wind(3);       d = wind(4);  % in the complex plane
% Increase 4 by 3 grid numbers for greater resolution
p = 400;  q = 300;               % p x q  grid of points

% p = 100;  q = 75;              % for rough drafts

one = ones(q,p);

x = a : (b-a)/(p-1) : b;

y = c : (d-c)/(q-1) : d;

[x,y] = meshgrid(x,y);           % grid of complex points
z = x + i*y;

clear x y

for n = 1:count
   s1 = z - r1;    s2 = z - r2;    s3 = z - r3;
   % value  = s1.*s2.*s3;
   % deriv = s2.*s3 + s1.*s3 + s1.*s2;
   % z = z - value./deriv;  	% Newton's iteration
   z = z - s1.*s2.*s3./(s2.*s3 + s1.*s3 + s1.*s2);
   n
end

clear value deriv s1 s2 s3

w1 = (abs(z - r1) < 0.01);    % Determine which root
w2 = (abs(z - r2) < 0.01);		% the initial guess z
w3 = (abs(z - r3) < 0.01);		% yields upon iteration 

w = 1 + w1 + 2*w2 + 3*w3;		% convert roots to colors

image(w);
colormap(newtmap);
axis off