xx=-4.0:0.01:4.0;
yy=xx;
[x,y]=meshgrid(xx,yy);
 
N=100;
z=x+i*y;
 
    for k=1:N
       z=z-(z.*z.*z-1.0)./(3.0.*z.*z);
 
% other functions:
%z=z-(z.*z.*z.*z-1.0)./(10.0);
%z=z-(z.*z-1.0)./(2.0.*z);
%z=z-(z.^4-z.^2+3.0)./(4.*z.^3-2.*z);
% z=z+(exp(-z)-1.0)./(exp(-z));
%z=z-(exp(-z)-z)./(1.0-exp(-z));
 
%z=z-(z.*z.*z.*z.*z.*z.*z.*z-1.0)./(8.0.*z.*z.*z.*z.*z.*z.*z);
%z=cos(z);
    end
 
imagesc(angle(z))
