function newtfractal(cs,xl,yl,n,nmax,tol)
% NEWTFRACTAL
% Function to demonstrate chaos in Newton's method
% Produces a colormap graph of how many iterations are required for
% Newton's method to converge to a root, 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, the number of Newton iterations are counted, up to a given maximum,
% and this is plotted according to color (using the defaul colormap, blue =
% immediate convergence, red = no convergence)
%
% syntax: newtfractal(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) maximum number of iterations -- default=200
%         tol    = (optional) error tolerance for Newton's method (ie
%                  iterate until f(z) < tol -- default=10^(-8)
%
% examples: newtfractal([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([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
%           taking anything above 50 iterations as non-convergent

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

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

deg = length(cs)-1;
dc = cs(1:deg).*(deg:-1:1);
q = nmax*ones(size(z));
fz = polyval(cs,z);
q(abs(fz)<tol) = 0;

for n=1:nmax
    z = z - fz./polyval(dc,z);
    fz = polyval(cs,z);
    q((abs(fz)<tol)&(q==nmax)) = n;
end


pcolor(x,y,q);
shading flat;
xlabel('Re(z_0)'), ylabel('Im(z_0)')
title('Number of iterations for Newton''s method to converge')