%
% qnewt.m   Steve Cox	cox@rice.edu
%
% color the newton basins for a quartic in a complex rectangle
%
% usage		qnewt(q,xt,yt,maxiter)
%
% where		q is the 5-vector of quartic coefficients
%		xt is the 3-vector, xlo xinc xhi
%		yt is the 3-vector, ylo yinc yhi
%		maxiter is the max number of iterations allowed
%
% example
%
%    qnewt([1 0 -0.84 0 -0.16],[.45 .0001 .55],[-.05 .0001 .05],20)
%

function qnewt(q,xt,yt,maxiter)

x = xt(1):xt(2):xt(3);
y = yt(1):yt(2):yt(3);
[X,Y] = meshgrid(x,y);
Z = X+i*Y;			% the complex matrix of newton seeds

dq = polyprime(q);		% the derivative of q

for k=1:maxiter,		% iterate every seed until maxiter

    Z = Z - polyval(q,Z)./polyval(dq,Z);

end

R = roots(q);				% roots of q

[r,j] = find(abs(Z-R(1))<0.1);		% find and plot the yellow seeds
plot(x(j),y(r),'y.','markersize',1)
hold on

[r,j] = find(abs(Z-R(2))<0.1);		% find and plot the red seeds
plot(x(j),y(r),'r.','markersize',1)

[r,j] = find(abs(Z-R(3))<0.1);		% find and plot the green seeds
plot(x(j),y(r),'g.','markersize',1)

[r,j] = find(abs(Z-R(4))<0.1);		% find and plot the blue seeds
plot(x(j),y(r),'b.','markersize',1)

hold off
axis off

qlab = '';		% initialize qlab, the string translation of q

for k=1:4,		% step through 1st 4 coefficients, checking
                        % for special cases (0,1,-1,i,-i)

    if q(k) ~= 0
       if isreal(q(k))		% real coefficient
          if q(k) == 1
             qlab = strcat(qlab,'z^',num2str(5-k),'+');
          elseif q(k) == -1
             qlab = strcat(qlab,'-z^',num2str(5-k),'+');
          else
             qlab = strcat(qlab,num2str(q(k)),'z^',num2str(5-k),'+');
          end
       else
          if real(q(k))==0 	% imaginary coefficient
             if imag(q(k)) == 1
                qlab = strcat(qlab,'iz^',num2str(5-k),'+');
             elseif imag(q(k)) == -1
                qlab = strcat(qlab,'-iz^',num2str(5-k),'+');
             else
                qlab = strcat(qlab,num2str(imag(q(k))),'iz^',num2str(5-k),'+');
             end
          else		% complex coefficient
             qlab = strcat(qlab,'(',num2str(q(k)),')z^',num2str(5-k),'+');
          end
       end
    end

end

if q(5) ~= 0		% take care of last coefficient
   if isreal(q(5))
      qlab = strcat(qlab,num2str(q(5)));
   else
      if real(q(5))==0
         qlab = strcat(qlab,num2str(imag(q(5))),'i');
      else
         qlab = strcat(qlab,num2str(q(5)));
      end
   end
end

qlab = strrep(qlab,'+-','-');		% replace +- with -
qlab = strrep(qlab,'-1i','-i');         % replace -1i with -i
qlab = strrep(qlab,'+1i','+i');         % replace +1i with +i

title(qlab,'fontsize',16)		% lay it on

% here is the replacememnt for polyder. One does not have to
% write a separate function, i.e., it is OK to just do this inline
% inside of qnewt

function dq = polyprime(q)		% differentiate the quartic
dq = [4 3 2 1].*q(1:4);