627.sagews

##example of Divided Differences
#(I should really make this into a method .. perhaps an exercise for my students :)) 

#points to interpolate
xPoints = [4,0,1,3]
yPoints = [-1, 3, 2, 3]

arr = [yPoints]

n = len(xPoints)  #number of points

#create divided differences table
for i in range(1,n):
    arr.append([])
    for j in range(n-i):
        calc = float((arr[i-1][j+1] - arr[i-1][j])/(xPoints[j+i]-xPoints[j]))
        arr[i].append(calc)

#coefficients
c = [arr[i][0] for i in range(n)]

x = var('x')

#create an array of the terms for interpolation error polynomial
terms = [c[0]];
for i in range(1,n):
    f = c[i]
    for j in range(0,i):
        f = f*(x-xPoints[j])
    terms.append(f)
   
#sum the terms together and place in a fuction
func = sum(terms)

points = zip(xPoints,yPoints)  #points to interpolate

p1 = list_plot(points, rgbcolor = (1,0,0))  #plot points

p = plot(func,min(xPoints)-1.0,max(xPoints)+1.0) #plot polynomial 
(p + p1).show(xmin=min(xPoints), xmax = max(xPoints), ymax=max(yPoints),ymin=min(yPoints))