%
% Draw an approximate boundary between basins of attraction
% in Newton iteration for f(x) = x^3 - 1
%
% RMC November 2010
%
f = @(x) x.^3-1;
df = @(x) 3*x.^2;
Newt = @(x) x - f(x)./df(x);
% The following takes quite a bit of memory 
% but makes a nice, lacy picture
lots = 1001;
x = linspace( -2, 2, lots );
y = x';
z = ones(lots,1)*x + 1i*y*ones(1,lots);
% Simple initial guess
r = z;
% Thirty iterations is plenty
for j=1:30,
  r = Newt(r);
end;
% The only points left are the roots or the boundary.
contour(x,y,r,'k');
axis('square');
set(gca,'fontsize',16);