function newtonfraktal(k)
%  Computes and plots the number of iterations in the complex 
%  newton iteration when solving the complex equation z^k=1.
%  It uses the subfunction  newtonc.m 
%  Written by Anders Logg 1999, modified by /stig 

figure(k)

x     = -1:0.01:1;
[X,Y] = meshgrid(x,x);
Z     = newtonc(X,Y,k);

surf(X,Y,Z)
view(2)
shading interp
axis square
title(['z^', num2str(k), ' = 1'])

%--subfunction----------------------------------------------------------

function n = newtonc(zr,zi,k)
% Computes the number of iterations in newton's method for solving the 
% complex equation z^k=1
% Written by Anders Logg 1999, modified by /stig

TOL_dz = 0.001;
TOL_n  = 20;  

n  = ones(size(zr));
dz = 1;

z = zr + i*zi;

while ( max(max(abs(dz))) > TOL_dz & max(max(n)) < TOL_n )
   dz = - (z.^k-1)./(k*z.^(k-1)+eps);  % z^k=1
   z  = z + dz;
   
   index = find(abs(dz) > TOL_dz);
   n(index) = n(index) + 1;

   
end	

