There's a good explanation starting on this web page.
I recommend that you read the three sections on Bezier surfaces in "Unit 8".
The process you described in your question is roughly correct.
Let's call the 16 points:
$$P_{11}, P_{12}, P_{13}, P_{14}$$
$$P_{21}, P_{22}, P_{23}, P_{24}$$
$$P_{31}, P_{32}, P_{33}, P_{34}$$
$$P_{41}, P_{42}, P_{43}, P_{44}$$
And let $\phi_1$, $\phi_2$, $\phi_3$, $\phi_4$ be the Bernstein polynomials of degree 3.
Suppose we want to evaluate a point on the surface with parameter values $(s,t) = (0.5, 0.2)$.
The first thing we do is use the 16 points four at-a-time to construct four curves. We start with the points of the first row: $P_{11}$, $P_{12}$, $P_{13}$, $P_{14}$. We can use them to construct a Bezier curve
$$
C_1(s) = \phi_1(s)P_{11} + \phi_2(s)P_{12} + \phi_3(s)P_{13} + \phi_4(s)P_{14} =
\sum_{j=1}^4 \phi_j(s)P_{1,j}
$$
Then, on this curve, we can calculate the point $Q_1 = C_1(0.5)$ at value $s=0.5$.
$$
Q_1 = C_1(0.5) = \sum_{j=1}^4 \phi_j(0.5)P_{1,j}
$$
You can calculate $Q_1$ by the normal de Casteljau algrorithm, or any other method.
We do the same things with the other three rows of points. For each row, we construct a Bezier curve, and then we calculate the point at parameter value $s=0.5$. We get points
$$\text{second row:}\quad Q_2 = C_2(0.5) = \sum_{j=1}^4 \phi_j(0.5)P_{2,j}$$
$$\text{third row:}\quad Q_3 = C_3(0.5) = \sum_{j=1}^4 \phi_j(0.5)P_{3,j}$$
$$\text{fourth row:}\quad Q_4 = C_4(0.5) = \sum_{j=1}^4 \phi_j(0.5)P_{4,j}$$
We can write all this more tidily as:
$$
Q_i = C_i(0.5) = \sum_{j=1}^4 \phi_j(0.5)P_{ij} \quad (i=1,2,3,4)
$$
Now we just do the Bezier thing in the other direction -- we use the points $Q_1$, $Q_2$, $Q_3$, $Q_4$ to construct a Bezier curve Q(t) given by
$$
Q(t) = \phi_1(t)Q_1 + \phi_2(t)Q_2 + \phi_3(t)Q_3 + \phi_4(t)Q_4 =
\sum_{i=1}^4 \phi_i(t)Q_i
$$
Then our final surface point is the point at parameter value $t=0.2$ on the curve $Q(t)$:
$$
X(0.5, 0.2) = Q(0.2) = \sum_{i=1}^4 \phi_i(0.2)Q_i
$$
Again, you can calculate this point using the de Casteljau algorithm or any other method.
To understand how all this relates to the surface equation, let's substitute the expressions for the $Q_i$ into our last equation. We get
$$
X(0.5, 0.2) = \sum_{i=1}^4 \phi_i(0.2)Q_i =
\sum_{i=1}^4 \phi_i(0.2)\left(\sum_{j=1}^4 \phi_j(0.5)P_{ij}\right)
$$
Rearranging this gives:
$$
X(0.5, 0.2) = \sum_{i=1}^4 \sum_{j=1}^4 \phi_j(0.5) \phi_i(0.2) P_{ij}
$$
You should recognize the right-hand side as the formula for a Bezier surface.
Here's a picture:
The blue curves are $C_1$, $C_2$, $C_3$, $C_4$, the red dots are $Q_1$, $Q_2$, $Q_3$, $Q_4$, and the red curve is $Q$.
Best Answer
Expand the parametric equation of the splitting curve. This is obtained by expliciting the Coons patch equation from the Bezier sides and freezing one of the parameters to $\frac12$.
Doing so, you will obtain a polynomial curve (if I am right, of degree one more than the initial Beziers). Then the X, Y and Z polynomials can be transformed to the Bernstein basis (by means of a linear transform), giving you the coordinates of the required control points.