%script file to plot newton iteration functions, using a 
%version of the escape time algorithm, coloring according to
%which root a value converges to and how quickly it does so

%create real and imaginary value arrays
[X,Y]=meshgrid([-2:0.01:2]);
[m,n]=size(X);
maxiters=30;
T=64*ones(m,n);
%parameters for the function C(z)=z^3+a*z+b
a=1.5;
b=0;
p=[1 0 a b];
%numerical values for the roots
rts=roots(p);
for j=1:m
    fprintf ('on iterate %d\n',j);
    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);
            z=y;
            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