%Stephen Muir's script for Complex Newton's Method
%  using f(z)=z^3+az+b (a, b are changable parameters)

%Create real and imaginary value arrays  These will be the initial points
%for Newton's Method- That is, we're dividing up the complex plane into
% an discrete array of points.
[X,Y]=meshgrid([-2:0.01:2]);
[m,n]=size(X);

maxiters=30;
% The matrix T will hold the color for each X+iY
T=64*ones(m,n);  %"64" is the default coloring for each point

%parameters for the function C(z)=z^3+a*z+b
a=0;
b=-1;
p=[1 0 a b];  %The coefficient vector for C(z)

%numerical values for the roots:  We use Matlab's built-in
% routine so we can see if Newton is converging to a root.
rts=roots(p);  %Matlab's root finding algorithm

for j=1:m
    fprintf ('Computing row %d out of %d\n',j,m);
    for k=1:n
        z=X(j,k)+i*Y(j,k);
        for t=1:maxiters
            y=z-(z^3+a*z+b)/(3*z^2+a);  %Newton's Method
            z=y;
            %Is our point close to a root?  If so, color its starting pt
            if abs(z-rts(1))<.001
                T(j,k)=1+t*(23/maxiters);
                break;
            end
            if abs(z-rts(2))<.001
                T(j,k)=30+t*(10/maxiters);
                break;
            end
            if abs(z-rts(3))<.001
                T(j,k)=45+t*(24/maxiters);
                break;
            end
        end
    end
end
image(T)
colorbar