import numpy as np
import pylab as pl
from matplotlib import cm

def iter_point(c):
    z=c
    tol=1.0e-2
    for i in range(1,1000):
        err=-(z**3-1)/(3*z**2)
        if abs(err)             break
        z=z-(z**3-1)/(3*z**2)
    return i

def draw_newton(cx,cy,d):
    x0,x1,y0,y1=cx-d,cx+d,cy-d,cy+d
    y,x=np.ogrid[y0:y1:200j,x0:x1:200j]
    c=x+y*1j
    newton=np.frompyfunc(iter_point,1,1)(c).astype(np.float)
    pl.imshow(newton,cmap=cm.Blues_r,extent=[x0,x1,y0,y1])
    pl.gca().set_axis_off()
    return

if __name__=='__main__':
    x,y=0.0,0.0
    draw_newton(x,y,0.5)
    pl.show()