By "column major convention," I assume you mean "The things I'm transforming are represented by $4 \times 1$ vectors, typically with a "1" in the last entry. That's certainly consistent with the second matrix you wrote, where you've placed the "displacement" $\Delta$ in the last column. Your entries in that second matrix follow a naming convention that's pretty horrible -- it's bound to lead to confusion.
Anyhow, the matrix product that you've computed is correct (that is, you did the multiplication properly). The result is something that first translates the origin to location $\Delta$ and the three standard basis vectors to the vectors you've called $\vec{x}$, $\vec{y}$, and $\vec{z}$, respectively, and having done so, then rotates the result in the $(2,3)$-plane of space (i.e., the plane in which the second and third coordinates vary, and the first is zero. Normally, I'd call this the $yz$-plane, but you've used up the names $y$ and $z$.) The rotation moves axis 2 towards axis 3 by angle $\theta$.
I don't know if that's what you want or not, but it's what you've got. Perhaps if you described the goal more clearly I could give a clearer answer.
If you begin with a horizontal object whose initial dimensions are $L_x$ and $L_y$ (before scaling), then the offset between blue and red text should be
$$
\begin{align}
&\Delta x={K_x-1\over2}L_x(1-\cos\alpha)+{K_y-1\over2}L_y\sin\alpha,\cr
&\Delta y=-{K_y-1\over2}L_y(1-\cos\alpha)+{K_x-1\over2}L_x\sin\alpha,\cr
\end{align}
$$
where $\alpha$ is the rotation angle with respect to the horizontal and $K_x$, $K_y$ are the scaling factors.
The sign of $\Delta x$ and $\Delta y$ depends on how you measure $\alpha$, so you may need to find the correct sign by trial and error. Hope that helps.
EDIT
To show how those formulas can be derived, consider the black rectangle in the diagram below. Scaling it (with respect to up-left vertex $P$) and then rotating around the center $O'$ of the scaled rectangle yields a final brown rectangle, whose up-left vertex is $P'$.
If the black rectangle is first rotated around its center $O$ and then scaled with respect to $P''$, then we get instead the blue rectangle below. We are interested in finding the "offset" vector $\vec{P'P''}$.
If we denote by $R$ the rotation operator
($R(x,y)=(\cos\alpha x-\sin\alpha y, \cos\alpha y+\sin\alpha x)$, if positive angles represent a counterclockwise rotation) and by $S$ the scaling operator
($S(x,y)=(K_x x,K_y y)$), then we have:
$$
\vec{P'P''}=\vec{P'O'}+\vec{O'O}+\vec{OP''}=
R\vec{PO'}+(S-1)\vec{OP}+R\vec{OP}=
-R(S\vec{OP})+(S-1)\vec{OP}+R\vec{OP},
$$
that is:
$$
\vec{P'P''}=
(S-1)\vec{OP}-R(S-1)\vec{OP}.
$$
Inserting here $\vec{OP}=(-L_x/2,L_y/2)$ one gets the desired result.
In comparing with my old formulas above, I see that
$\vec{P'P''}=(-\Delta x,-\Delta y)$.
In addition, I used there the opposite convention for the sign of the angle, but I hope this explanation is clear enough.
Best Answer
If you already implemented translation, then you can translate $P$, the center of your scaling, to $\underline{0}$, the center of your reference system. (The origin of the axis).
Then you apply a scaling centered in the origin: a scaling by a factor $k$ is represented by the matrix: $$\left[ \begin{matrix} k & 0 & 0 & 0 \\ 0 & k & 0 & 0 \\ 0 & 0 & k & 0 \\ 0 & 0 & 0 & 1 \end{matrix} \right] $$
And then you apply the opposite translation to what you did before, bringing $\underline{0}$ back to $P$.
If you want to use a single matrix transformation, then all you have to do is multiply the appropriate matrixes "by hand":
$$ \left[ \begin{matrix} 1 & 0 & 0 & P_1 \\ 0 & 1 & 0 & P_2 \\ 0 & 0 & 1 & P_3 \\ 0 & 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} k & 0 & 0 & 0 \\ 0 & k & 0 & 0 \\ 0 & 0 & k & 0 \\ 0 & 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} 1 & 0 & 0 & -P_1 \\ 0 & 1 & 0 & -P_2 \\ 0 & 0 & 1 & -P_3 \\ 0 & 0 & 0 & 1 \end{matrix} \right] = \left[ \begin{matrix} k & 0 & 0 & (1-k)P_1 \\ 0 & k & 0 & (1-k)P_2 \\ 0 & 0 & k & (1-k)P_3 \\ 0 & 0 & 0 & 1 \end{matrix} \right] $$