# given data points in an x-vector and y-vector below.
# this will compute the Newton finite difference 
# polynomial interpolation of those points
xx = [0,2,4,6,8,10]
yy = [10,9,4,-2,6,7]
var('x')

N = len(xx)
# first we need to make a finite difference table
#  We'll make a list of lists

Dyy = [yy]

# practice
#Dyy += [[yy[j+1]-yy[j] for j in range(N-1)]]
# practice #2 - exactly the same
#Dyy += [[Dyy[0][j+1]-Dyy[0][j] for j in range(N-1)]]
# Above we have added Dyy[1] to Dyy

# Now we'll do all of them
for k in range(1,N):
    m = len(Dyy)-1
    Dyy += [[Dyy[m][j+1]-Dyy[m][j] for j in range(N-1-m)]]

show(Dyy)
h = xx[1]-xx[0]

Newt(x) = sum([Dyy[j][0]/factorial(j)/h^j*prod([x-xx[k] for k in range(j)]) for j in range(N)])
show(Newt(x))
show(expand(Newt(x)))

pts = [[xx[j],yy[j]] for j in range(N)]
p1 = scatter_plot(pts,marker='o')
p2 = plot(Newt(x),(0,10),ymin=-2,ymax=10)
show(p1+p2)