reset()
xs = [-1, 0, 2]
ys = [2,-1,1]
m = [1, 1, 1]
D = zero_matrix(QQ,7)
# D.is_immutable()
DD = copy(D)
DD[0,0] = xs[0]
DD[1,0] = xs[0]
DD[2,0] = xs[1]
DD[3,0] = xs[1]
DD[4,0] = xs[2]
DD[5,0] = xs[2]
DD[1,2] = m[0]
DD[0,1] = ys[0]
DD[1,1] = ys[0]
DD[3,2] = m[1]
DD[2,1] = ys[1]
DD[3,1] = ys[1]
DD[5,2] = m[2]
DD[4,1] = ys[2]
DD[5,1] = ys[2]
for k in range(5):
if DD[k+1,0]-DD[k,0] != 0:
DD[k+1,2] = (DD[k+1,1]-DD[k,1])/(DD[k+1,0]-DD[k,0])
end
end
for k in range(1,5):
DD[k+1,3] = (DD[k+1,2]-DD[k,2])/(DD[k+1,0]-DD[k-1,0])
end
for k in range(2,5):
DD[k+1,4] = (DD[k+1,3]-DD[k,3])/(DD[k+1,0]-DD[k-2,0])
end
for k in range(3,5):
DD[k+1,5] = (DD[k+1,4]-DD[k,4])/(DD[k+1,0]-DD[k-3,0])
end
for k in range(4,5):
DD[k+1,6] = (DD[k+1,5]-DD[k,5])/(DD[k+1,0]-DD[k-4,0])
end
show(DD)
c = [DD[0,1], DD[1,2], DD[2,3], DD[3,4], DD[4,5], DD[5,6]]
show(c)
p = c[0] + (x-xs[0]) *( c[1] + (x-xs[0]) * (c[2] +(x-xs[1])*(c[3] + (x-xs[1])*(c[4] + (x-xs[2])*(c[5] ) ) )))
show(p)
G = plot(p,(x,-1,2))
G += point([(xs[0],ys[0]),(xs[1],ys[1]),(xs[2],ys[2])],size=30, color='red')
G += plot(m[0]*(x-xs[0])+ys[0],(x,xs[0]-0.2,xs[0]+0.2),color='green')
G += plot(m[1]*(x-xs[1])+ys[1],(x,xs[1]-0.2,xs[1]+0.2),color='green')
G += plot(m[2]*(x-xs[2])+ys[2],(x,xs[2]-0.2,xs[2]+0.2),color='green')
G.show(aspect_ratio=1,figsize=(6,8))
