interpolating using bezier splines

We will often want to use splines to model other data.  There are two modeling, interpolation and approximation. 

Wikipedia states: interpolation is a method of constructing new data points within the range of a discrete set of known data points.

A cubic Bezier spline segment, as a cubic polynomial can perfectly interpolate four points.

Finding the interpolating cubic polynomial for a set of four data points is a special case of more general polynomial interpolation.

The most direct form of this polynomial interpolation is to directly plug into this Lagrange basis functions for y=y[0],y[1],y[2],y[3] at t=0,1/3,2/3,1.

In this way, y = y[0] + t^2 (9 y[0] - (45 y[1])/2 + 18 y[2] - (9 y[3])/2) +
 t (-((11 y[0])/2) + 9 y[1] - (9 y[2])/2 + y[3]) +
 t^3 (-((9 y[0])/2) + (27 y[1])/2 - (27 y[2])/2 + (9 y[3])/2) and similarly for x.

When we expand into standard polynomial form we get:
x[t_]:= a_x + t * b_x + t^2 * c_x + t^3 *d_x
y[t_]:= a_y + t * b_y + t^2 * c_y + t^3 *d_y

giving

a_y = y[0]
b_y = (-((11 y[0])/2) + 9 y[1] - (9 y[2])/2 + y[3])
c_y = (9 y[0] - (45 y[1])/2 + 18 y[2] - (9 y[3])/2)
d_y = (-((9 y[0])/2) + (27 y[1])/2 - (27 y[2])/2 + (9 y[3])/2)

When we put this into the form of Bezier control points:

for f[0]={x[0],y[0]}, {x[1],y[1]}, {x[2],y[2]}, {x[3],y[3]}

p_0 = f[0]
p_1 = 1/6 (-5 f[0]+18 f[1]-9 f[2]+2 f[3])
p_2 = 1/6 (2 f[0]-9 f[1]+18 f[2]-5 f[3])
p_3 = f[3]