function newtonfractal(iter,npts)
% NEWTONFRACTAL plot out the fractals of Newton's method
%    This program concerns about Newton's method for finding the 
%    roots of the cubic polynomial z^3-1. Newton's method depends             
%    very sensitively on the initial guess. And the dependence is       
%    actually of fractal shape for many polynomials. The program returns           
%    a color picture where each color represents a particular root           
%    of the cubic polynomial z^3-1. For instance, if you start from a            
%    point in yellow, you will end up with the root z=1.
%    try, e.g.,  newtonfractal(50,200)

xaxis=0;
yaxis=0;
d=1.5;
x=linspace(xaxis-d,xaxis+d,npts);
y=linspace(yaxis-d,yaxis+d,npts);
[xtrans,ytrans]=meshgrid(x,y);
ztrans=xtrans+i*ytrans;

% run Newton's method sufficiently many times to get convergence for
% all starting points.
for k=1:iter;

% change this line for different polynomials
ztrans=2/3.*ztrans+1/3./ztrans.^2;
end

% imaginary parts of these three roots are all different and can thus
% be used to differentiate each root.
% change this line for different polynomials
t=imag(ztrans);

% plot out
colormap prism(256)
pcolor(t);
shading flat;
axis('square','equal','off');