#Cubic splines

x = [0,2,4,6,8,10]
y = [10,9,4,-2,6,7]
n = len(x)-1
A = list(var('A%d' % i) for i in range(n)) # this makes n variables A[0]=A0,A1,..
B = list(var('B%d' % i) for i in range(n))
C = list(var('C%d' % i) for i in range(n))
D = list(var('D%d' % i) for i in range(n))
eqs = [y[j]==A[j]*x[j]^3+B[j]*x[j]^2+C[j]*x[j]+D[j] for j in range(n)]
eqs += [y[j+1]==A[j]*x[j+1]^3+B[j]*x[j+1]^2+C[j]*x[j+1]+D[j] for j in range(n)]
eqs += [3*A[j]*x[j+1]^2+2*B[j]*x[j+1]+C[j]==3*A[j+1]*x[j+1]^2+2*B[j+1]*x[j+1]+C[j+1] for j in range(n-1)]
eqs += [6*A[j]*x[j+1]+2*B[j]==6*A[j+1]*x[j+1]+2*B[j+1] for j in range(n-1)]
eqs += [6*A[0]*x[0]+2*B[0]==0,6*A[n-1]*x[n]+2*B[n-1]==0]
vs = A+B+C+D
answer = solve(eqs,vs)
show(answer)
AA = [answer[0][j].rhs() for j in range(n)]
BB = [answer[0][j].rhs() for j in range(n,2*n)]
CC = [answer[0][j].rhs() for j in range(2*n,3*n)]
DD = [answer[0][j].rhs() for j in range(3*n,4*n)]

pts = [[x[j],y[j]] for j in range(n+1)]
p1 = scatter_plot(pts, marker='o')
var('xx')
bubba = p1
pp=[]
for j in range(n):
    pp=plot(AA[j]*xx^3+BB[j]*xx^2+CC[j]*xx+DD[j],(x[j],x[j+1]))
    bubba = bubba +pp
show(bubba)