matrix([1,-1],[2,4],[3,6],[4,3],[5,6]);
pts:transpose(%);

n:length(pts[1]);
f[k]:=sum(a[k,j]*t^j,j,0,3);
create_list(create_list(a[k,j],j,0,3),k,1,n-1);
vars:flatten(%);

eq:[]$
for i:1 thru n-1 do
eq:append(eq,[subst(t=pts[1,i],f[i])=pts[2,i],
        subst(t=pts[1,i+1],f[i])=pts[2,i+1]]
        )$
for i:2 thru n-1 do
    eq:append(eq,[subst(t=pts[1,i],diff(f[i-1],t))=
        subst(t=pts[1,i],diff(f[i],t)),
        subst(t=pts[1,i],diff(f[i-1],t,2))=
        subst(t=pts[1,i],diff(f[i],t,2))])$
eq:append(eq,[subst(t=pts[1,1],diff(f[1],t,2))=0,
        subst(t=pts[1,n],diff(f[n-1],t,2)=0)])$
sol:solve(eq,vars);

eq;

for i:1 thru n-1 do
    fa[i]:substitute(sol,f[i])$
pf:create_list([parametric,t,fa[i],[t,pts[1,i],pts[1,i+1]]],i,1,n-1)$

plot2d(pf)$