Scaling a rotated child with parent-child transforms

Question

In unity, if you parent something to another object, then give the child some rotation, and scale the parent, the scale applies to some arbitrary axis to the child, here's some pictures to demonstrate:

Now I scale the Z axis of the parent object and this happens:

My from scratch game does not do this. Rather it applies the parent scale along the same local axis on the child, so you get this:

I don't know how to achieve Unity's scaling behaviour though. Right now transformations are combinations of standard transform matrices that are just applied child to parent. Like this:

Matrixf child = Matrixf::Translate(childTrans) * Matrixf::Rotate(childRotation) * Matrixf::Scale(childScale);
Matrixf parent = Matrixf::Translate(parentTrans) * Matrixf::Rotate(parentRotation) * Matrixf::Scale(parentScale);

Matrixf worldChildTransform = parent * child;

Transform matrices defined as such:

inline static Matrix Translate(Vec3<T> translate)
{
    Matrix mat;
    mat.m[0][0] = T(1.0); mat.m[0][1] = T(0.0); mat.m[0][2] = T(0.0); mat.m[0][3] = translate.x;
    mat.m[1][0] = T(0.0); mat.m[1][1] = T(1.0); mat.m[1][2] = T(0.0); mat.m[1][3] = translate.y;
    mat.m[2][0] = T(0.0); mat.m[2][1] = T(0.0); mat.m[2][2] = T(1.0); mat.m[2][3] = translate.z;
    mat.m[3][0] = T(0.0); mat.m[3][1] = T(0.0); mat.m[3][2] = T(0.0); mat.m[3][3] = T(1.0);
    return mat;
}

inline static Matrix Rotate(Vec3<T> rotation)
{
    // This is a body 3-2-1 (z, then y, then x) rotation
    const T cx = cos(rotation.x);
    const T sx = sin(rotation.x);
    const T cy = cos(rotation.y);
    const T sy = sin(rotation.y);
    const T cz = cos(rotation.z);
    const T sz = sin(rotation.z);

    Matrix res;
    res.m[0][0] = cy*cz;    res.m[0][1] = -cx*sz + sx*sy*cz;    res.m[0][2] =  sx*sz + cx*sy*cz;    res.m[0][3] = T(0.0);
    res.m[1][0] = cy*sz;    res.m[1][1] =  cx*cz + sx*sy*sz;    res.m[1][2] = -sx*cz + cx*sy*sz;    res.m[1][3] = T(0.0);
    res.m[2][0] = -sy;      res.m[2][1] = sx*cy;                res.m[2][2] = cx*cy;                res.m[2][3] = T(0.0);
    res.m[3][0] = T(0.0);   res.m[3][1] = T(0.0);               res.m[3][2] = T(0.0);               res.m[3][3] = T(1.0);
    return res;
}

inline static Matrix Scale(Vec3<T> scale)
{
    Matrix mat;
    mat.m[0][0] = scale.x;  mat.m[0][1] = T(0.0);   mat.m[0][2] = T(0.0);   mat.m[0][3] = T(0.0);
    mat.m[1][0] = T(0.0);   mat.m[1][1] = scale.y;  mat.m[1][2] = T(0.0);   mat.m[1][3] = T(0.0);
    mat.m[2][0] = T(0.0);   mat.m[2][1] = T(0.0);   mat.m[2][2] = scale.z;  mat.m[2][3] = T(0.0);
    mat.m[3][0] = T(0.0);   mat.m[3][1] = T(0.0);   mat.m[3][2] = T(0.0);   mat.m[3][3] = T(1.0);
    return mat;
}

EDIT: Some more information about the matrices in this specific case:

Given this configuration

Parent position = (0, 0, 3)
Parent scale = (0.5, 0.5, 0.2)
Parent rotation = (0, 0, 0) (euler angles, radians)

Child local position = (0, 0, 4.2)
Child local scale = (1, 1, 1)
Child local rotation = (0, 0.9, 0) (euler angles, radians)

Then we have the following matrices

Child (combined translation, rotation and scale)

{0.621609986, 0.000000000, 0.783326924, 0.000000000}
{0.000000000, 1.00000000, 0.000000000, 0.000000000}
{-0.783326924, 0.000000000, 0.621609986, 4.19999981}
{0.000000000, 0.000000000, 0.000000000, 1.00000000}

Parent (combined again)

{0.500000000, 0.000000000, 0.000000000, 0.000000000}
{0.000000000, 0.500000000, 0.000000000, 0.000000000}
{0.000000000, 0.000000000, 0.200000003, 3.00000000}
{0.000000000, 0.000000000, 0.000000000, 1.00000000}

Then parent * child gives this:

{0.310804993, 0.000000000, 0.391663462, 0.000000000}
{0.000000000, 0.500000000, 0.000000000, 0.000000000}
{-0.156665385, 0.000000000, 0.124321997, 3.83999991}
{0.000000000, 0.000000000, 0.000000000, 1.00000000}

Answer 1

After some experiments, I figured this out, and as was noticed in the comments to my question it was not to do with my matrix multiplication order, but rather something else entirely.

After calculating the world transform of my object, I decomposed it into a vec3 translation, vec3 scale and quat rotation. This operation stripped out what is effectively a shear operation on the child cube. Now I just save the world matrix for use in rendering later, and the problem is solved!