How do I calculate the values of the control points for an uniform cubic B spline surface?

How do I calculate the values of the control points for an uniform cubic B
spline surface?

I want to interpolate the following 3 scattered data points:
(80.9,58.5,48.0),(35.0,89.6,82.3),(74.7,17.4,85.9) by an uniform cubic B
spline surface on the following control lattice, $ \phi $:
In this figure the data points (blue) have been normalized to [0,1]: $ x_n
= (x-xmin)/(xmax-xmin) $
The B spline interpolator becomes:
$ z(x,y) = \Sigma_{i=-1}^2\Sigma_{j=-1}^2B_{i+1}(x_n)B_{j+1}(y_n)\phi_{i,j} $
where
$ B_0 (t)=(1-t^3)/6 $
$ B_1 (t)=(3t^3-6t^2+4)/6 $
$ B_2 (t)=(-3t^3+3t^2+3t+1)/6 $
$ B_3 (t)=t^3/6 $
How do I calculate the values of the control points $ \phi_{i,j} $?
Introducing:
$ \vec{z} = [z_{0,0},z_{0,1},...,z_{3,3}]^T $,
$ A_i = [B_0(x_i)B_0(y_i), B_0(x_i)B_1(y_i), ..., B_3(x_i)B_3(y_i)] $,
$ A=[A_1^T, A_2^T, A_3^T]^T $,
$ \vec{\phi} = [\phi_{0,0},\phi_{0,1},...,\phi_{3,3}]^T $
The interpolation can be written as a matrix equation:
$ \vec{z} = A\vec{\phi} $
I can then find a Least Square Error solution to this equation.
Doing this in Matlab and sampling the resulting interpolator gives the
following result:
The interpolator goes exactly trough the data points (black) but the
extrapolation in the lower left corner seems bad.
What is causing this problem?
I would guess the matrix equation is underdetermined (3 equations and 16
unknowns). Is there some easy way I can regularize it so that the
interpolator gives me a more planar solution?