def find_roots(method,a,b,segments=10):
    seg = np.linspace(a,b,segments+1);
    for k in range(segments):
        ak, bk = seg[k:k+2];
        #print "searching for roots in",[ak,bk]
        x = (ak+bk)/2;
        count = 0;
        while ak<=x<=bk and count < 50:
            count += 1;
            xold, x = x, method(x);
            #print x
            if count==2 and abs(x-xold)>1e-1*(bk-ak): break;
            if abs(x-xold)<1e-8:
                y,_,_ = f(x)
                print "found root x=%.15f with f(x)=%.8e in %d iterations"%(x,y,count);
                break;

find roots with Newton step
found root x=3.000000007315551 with f(x)=-3.77475828e-15 in 23 iterations
found root x=10.999999991701889 with f(x)=-3.33066907e-16 in 23 iterations
find roots with Halley step
found root x=3.000000004913715 with f(x)=-1.66533454e-15 in 15 iterations
found root x=10.999999999234854 with f(x)=0.00000000e+00 in 16 iterations
find roots with Newton plus double Newton step
found root x=2.999999999980970 with f(x)=0.00000000e+00 in 4 iterations
found root x=10.999999999997232 with f(x)=0.00000000e+00 in 3 iterations