% script file Ex8.m
% compute roots of polynomial z^3-1 with
% different initial values
% last changes: 17.11.01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INITIALIZATION

% prepare matrix of initial values
nx = 1001; ny = 1001;     % number of grid points
ax = -1;  bx = 1;
ay = -1;  by = 1;
x  = linspace(ax,bx,nx); % real value
y  = linspace(ay,by,ny); % imaginary value

[X,Y] = meshgrid(x,y);
Z     = X+i*Y;

itmax = 60;   % maximal number if iterations
toler = 1e-5; % decision tolerance
diver = 100;   % "divergence"

% analytically known roots:
w = [1 , exp(2*pi*i*1/3) , exp(2*pi*i*2/3)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ITERATION

% Simply iterate each component itmax times. More smart would be to 
% put a stopping criterion as "many" components converge "fast". 
% however, this would be of importance if the iteration number
% which is needed would be high. For the moment convergence check
% in costs is camparable with costs for more iterations ...

Z2 = zeros(size(Z));

for k = 1:itmax
   Z1 = Z - (Z.^3-1)./(3*Z.^2);
   Z = Z1;  
end

% decide to which w the sequence has converged
Limit = (abs(Z1-w(1))<toler)*1 + (abs(Z1-w(2))<toler)*2 + (abs(Z1-w(3))<toler)*3 + ...
   (abs(Z1)>diver|isnan(Z1))*4;

% after this construction, matrix Limit has values
% 1 at points which have converged to w(1)
% 2                                   w(2)
% 3                                   w(3)
% 4                      diverged
% 0 at all other points (i.e., not yet converged).

% if we have points where Limit has zero values,
% itmax must be increased, as the convergence is not yet
% achieved

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOTTING

imagesc(x,y,flipud(Limit)) % that's bit faster than surf, but has 
                         % that funny MATLAB flipud problem...
xlabel('Re(z)')
ylabel('Im(z)')

% tell to use only 4 or 5 colours in plots

if ~isempty(find(Limit==0))
   warning('Convergence was not yet achieved. Increase parameter itmax!')
	MAP = [1 1 1; ... % white for points which have not yet converged
   	   1 0 0; ...  % red
      	0 1 0; ...  % green
      	0 0 1; ...  % blue
         0 0 0];     % black
else
   MAP = [1 0 0; ... % red
      	0 1 0; ...  % green
      	0 0 1; ...  % blue
         0 0 0];     % black
end
colormap(MAP)

% mark the stationary points
hold on
P=plot3(real(w),imag(w),5*ones(size(w)),'y.'); % factor 5 needed to move it above the 
                                               % previous plot
set(P,'Markersize',10)
view(2)
axis square

