How to calculate the control points of a Bézier curve?

Here is my approach:

pts = {{0, 0}, {2, 1}, {4, 3}, {6, 1}};
paras = FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts]), Total]] // N
mat = Outer[BernsteinBasis[3, #1, #2] &, Range[0, 3], paras] // Transpose;
ctrlpts = LinearSolve[mat, pts]
(* {{0., 0.}, {2.71043, -0.262717}, {3.94236, 6.32778}, {6., 1.}} *)

Graphics[{BezierCurve[ctrlpts], PointSize[Medium], Red, Point[pts]}]

Comparison with Interpolation[]

To show that the results of the Bézier curve interpolant and the built-in Interpolation[] are vastly different, I will use the following data:

pts1 = {{-1, 0}, {2, 1}, {4, 4}, {6, -3}};
mat1 = Outer[BernsteinBasis[3, #1, #2] &, Range[0, 3], 
             FoldList[Plus, 0.0, 
                      Normalize[(Norm /@ Differences[pts1]), Total]]] // Transpose;

f = Interpolation[pts1];
Show[Plot[f[x], {x, -1, 6}], 
     Graphics[{BezierCurve[LinearSolve[mat1, pts1]], 
               PointSize[Medium], Red, Point[pts1]}], 
     PlotRange -> {Automatic, {-3, 6}}]

There's no need to resort to chord length parametrization as was done in the other answer for this simple interpolation problem.

pts = {{0, 0}, {2, 1}, {4, 3}, {6, 1}};
{xt, yt} = Transpose[pts];
cp = LinearSolve[Outer[BernsteinBasis[3, #2, #1] &, Rescale[xt], Range[0, 3]], yt]
   {0, -7/6, 20/3, 1}

Graphics[{{Directive[AbsoluteThickness[2], ColorData[97, 1]], 
           BezierCurve[Transpose[{xt, cp}]]},
          {Directive[AbsolutePointSize[6], ColorData[97, 2]], Point[pts]}}]

An explanation for the following two results is left as an exercise for the interested reader:

if = Interpolation[pts];
Plot[if[x], {x, 0, 6}, 
     Epilog -> {Directive[AbsoluteThickness[2], ColorData[97, 2]],
                BezierCurve[Transpose[{xt, cp}]]}, 
     PlotStyle -> AbsoluteThickness[6]]

Plot[if[x] - cp.BernsteinBasis[3, Range[0, 3], Rescale[x, MinMax[xt]]],
     {x, 0, 6}, Evaluated -> True, PlotRange -> All]