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

N = len(xx)

L(x) = sum([yy[j]*prod( [(x-xx[k])/(xx[j]-xx[k]) for k in range(N) if j!=k]) for j in range(N)])
show(L(x))
show(expand(L(x)))
pts = [[xx[j],yy[j]] for j in range(N)]
p1 = scatter_plot(pts,marker='o')
p2 = plot(L(x),(-2,11),ymin=-25,ymax=25)
show(p1+p2)