function newtfractal_boa(cs,xl,yl,n,nmax,tol)
% NEWTFRACTAL
% Function to demonstrate chaos in Newton's method
% Produces a colormap graph of the basins of attraction for Newton's method,
% as a function of initial starting position z_0.  The function for which
% the roots are being found is a complex polynomial and z_0 is a point in
% the complex plane.  For each z_0, a fixed number of iterations of
% Newton's method is applied and the root to which z_0 converges is
% determined; this is then plotted according to color -- all z_0 that
% converge to one root are given one color, those that converge to another
% root are given a second color, and so on.  In the default colormap, red
% corresponds to points that do not converge at all (in the given number of
% iterations).
%
% syntax: newtfractal_boa(coeffs,[xl,xu],[yl,yu],n,Nmax,tol)
%
% inputs: coeffs = vector of polynomial coefficients
%         xl, xu = lower and upper bounds for the real part of z_0
%         yl, yu = lower and upper bounds for the imaginary part of z_0
%         n      = number of gridpoints to use in x and y
%         Nmax   = (optional) number of Newton iterations -- default=100
%         tol    = (optional) error tolerance for determining whether
%         Newton's method converged (ie convergence if z_Nmax is within
%         tol of a root) -- default=10^(-6)
%
% examples: newtfractal_boa([4,0,0,0,-1],[-1,1],[-2,2],151)
%           finds roots of z^4 - 1 = 0 (ie 4th roots of 1, which are +/- 1
%           and +/- i) with Re(z_0) in [-1,1] and Im(z_0) in [-2,2] with
%           151 points in each direction
%
%           newtfractal_boa([4,0,0,0,-1],[-0.25,0.25],[-0.25,0.25],201,50)
%           same as above but on a smaller region, with more points, and
%           using only 50 iterations

if (nargin<6)
    tol=1e-6;
    if (nargin<5)
        nmax=100;
    end
end
x=linspace(xl(1),xl(2),n);
y=linspace(yl(1),yl(2),n);

[u,v]=meshgrid(x,y);
z = u+i*v;

deg = length(cs)-1;
dc = cs(1:deg).*(deg:-1:1);

for n=1:nmax
    z = z - polyval(cs,z)./polyval(dc,z);
end

r=roots(cs);
q = deg * ones(size(z));
for k = 1:deg
   index = abs(z - r(k)) < tol;
   q(index) = k-1;
end

pcolor(x,y,q);
shading flat;
xlabel('Re(z_0)'), ylabel('Im(z_0)')
title('Basins of attraction for Newton''s method')