CubicSpline := proc(  x :: Vector,
                      y :: Vector,
                      t :: name
                     )
   local n, p, k, m, eq, a;

   # get the number of points
   n := LinearAlgebra:-Dimension(x);

   a := 'a';

   # define the spline polynomials
   for k from 1 to n-1 do
      p[k] := a[k,1] + a[k,2]*(t-x[k]) 
             + a[k,3]*(t-x[k])^2 + a[k,4]*(t-x[k])^3;
   end do;

   # now define 4*n-4 equations
   m := 0;  # equation counter

   # obvious equations
   for k from 1 to n-1 do
      m := m + 1; 
      eq[m] := subs(t=x[k],p[k]) = y[k];
      m := m + 1;
      eq[m] := subs(t=x[k+1],p[k]) = y[k+1];
   end do;

   # join the first derivative
   for k from 1 to n-2 do
      m := m + 1;
      eq[m] := subs(t=x[k+1],diff(p[k],t)) 
               = subs(t=x[k+1],diff(p[k+1],t));
   end do;

   # join the second derivative
   for k from 1 to n-2 do
      m := m + 1;
      eq[m] := subs(t=x[k+1],diff(p[k],t$2)) 
               = subs(t=x[k+1],diff(p[k+1],t$2));
   end do;

   m := m + 1;
   eq[m] := subs(t=x[1], diff(p[1],t$2) ) = 0;
   m := m + 1;
   eq[m] := subs(t=x[n], diff(p[n-1],t$2) ) = 0;

   # solve equations
   solve( {seq(eq[k],k=1..m)}, 
         {seq( seq(a[k,j],j=1..4), k=1..n-1) } );
   assign(%);

   # output
   return piecewise( seq(op([t<=x[k+1],p[k]]),k=1..n-1) );
end proc;