newton.py

from numpy.polynomial.polynomial import polyval
from scipy.optimize import newton
from random import random

def poly(coeffs):
    def polyfunc(x):
        return polyval(x, coeffs)
    return polyfunc

# Polynomial: coeffs[0]*x^0 + coeffs[1]*x^1 + coeffs[2]*x^2 + ... 
# coeffs are random numbers between 0 and 1
coeffs = [random() for i in range(41)]

print 'y = ' + ' + '.join(['%.3f*x^%d' % (coeff, i) for i, coeff in enumerate(coeffs)])
print 'x0 = 0'

print newton(poly(coeffs), 0)