w = 0.5;
n = 100;
divergence = 1e6;  % condition for divergence of iteration 
root_tol   = 1e-5;
it_tol     = 1e-8;

% x goes from -w to w, y goes from -iw to iw
x = linspace(-w, w, n);
y = ...
% create the section of complex plane with the above bounds
z = ...

% this function is just f(z) =  z^3 -1
f  = @(mat) mat.^3 -1;
fp = ...

newton_it = ...

% now lets perform this newton iteration on the complex plane until things
% stop changing
it_diff = 1;
while it_diff > it_tol                     % while not converged
    z_new = ...                            % perform newton iteration
    z_new(...) = NaN;                      % identify values greater than divergence and set to NaN for identification
    diff = z_new - z;                      % subtract new values from old values
    diff(isnan(diff)) = 0;                 % all values that diverged aren't changing anymore
    it_diff = norm(diff);                  % take norm of the diff matrix to determine stopping condition
    
    z = z_new;                             % set z to the new value
end

% now we want to assign a unique number to each root that we found 

% look at all of the values in z as a single vector
roots = z(:);
% sort values for easy comparison
roots = sort(roots,'ComparisonMethod','real');
% compare each value to it's neighbour
roots_shift = circshift(roots,1);
% unique roots correspond to locations where neighbours change
unique_locations = roots - roots_shift;

% get our unique roots
unique_roots = roots((abs(unique_locations) > root_tol) | (isnan(unique_locations)));
colors = z(:,:);

% loop through and assign unique number to each type of root
for ii=1:length(unique_roots)
    colors(abs(z-unique_roots(ii)) <= root_tol) = ii;
end
colors(isnan(colors)) = length(unique_roots);

figure;
sub1=subplot(1,1,1); 
colormap jet;
image(x,x,colors, 'CDataMapping','scaled')
title('roots!');

