Fractal solution, Plotting complex plane in Matlab

The following is a solution set for a function f(z) =z^4+8iz^2-25. The roots of the function are
2-i, -2+i, -1+2i and 1-2i
The color shows which values will converge to the solution in the complex place


The complete matlab program is as below


% Using Newton Basin method to plot in complex plane
% The complex function is z^4+8iz^2-25
% There are 4 roots 2-i, -2+i, -1+2i and 1-2i of the above function
% It shows how starting with different points root converge

tol =.01;
a=0;
b=2;
c=0;
d=2;
m=1000;
n=1000;
x=linspace(a,b,m); % Setting up the x axis Coordinate System a to b
y=linspace(c,d,n); % Setting up the y axis Coordinate System c to d
z1=2-i; % root 1
z2=-2+i; % root 2
z3=-1+2i; % root 3
z4=1-2i; % root 4
lmax=20;
r=ones(n,m); % Setting a n*m matrix to be all 1
for j=1:n
for k=1:m
z=x(k)+i*y(n-j+1); % Generating the complex number, Note n-j+1
if z == 0;
z = tol;
endif

for l = 1:lmax % Newton iteration
zz = (3*z^4+8*i*z^2+25)/(4*z^3+16*i*z); % Newton's Formula for Iteration
if abs(z-zz) < tol
if abs(z-z1) < tol
r(j,k) =1; %Red
elseif abs(z-z2) r(j,k) =55; %Pink
elseif abs(z-z3) r(j,k) =33; %Light Blue
elseif abs(z-z4) < tol
r(j,k) = 42; %Deep Blue
endif
break
else
z=zz;
endif

endfor
endfor
endfor

colormap('hsv')
image(r)
axis square
axis off