1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:

clear all; 
format compact; 
z=[1 0 0 -1] 
zi=roots(z) 
nmax=20 
ab=.0001 
for a=-2:.1:2 
    for b=-2:.1:2 
        zn=complex(a,b); 
        for k=1:nmax 
               if abs(zn-zi(1))*abs(zn-zi(2))*abs(zn-zi(3)) > ab 
                  zn=zn-(zn^3-1)/(3*zn^2+eps) 
                else 
                    k=nmax; 
                    if zn-zi(1) < ab; 
                        plot(a,b,'r'); 
                    elseif zn-zi(2) < ab; 
                        plot(a,b,'g'); 
                    elseif zn-zi(3) < ab; 
                        plot(a,b,'b'); 
                    end 
                end 
            end 
        end 
    end