Justification for the post multiplication of Euler angle matrices.

The basic matrix representations of rotations about the three fixed axes of an orthonormal coordinate system are easy to derive.
For right-handed rotation of an angle a, and column vectors, we get the following:

Around x,    Around y,    Around z
| 1 0 0  |   | c 0 s  |   | c -s 0 |
| 0 c -s |   | 0 1 0  |   | s c 0  |
| 0 s c  |   | -s 0 c |   | 0 0 1  |

Where c = cosine(a), s = sin(a).

Note: An easy way to derive these matrices without having to go back to pictures is use 'c' for the direction cosine of a component along it's direction, and 's' or '-s' along the remaining direction depending on the order of the components in the string x,y,z,x,y,z.
So per example, x on y gives 's' (since we have x,y in the string) while y on x gives '-s' (since we have z in between: y,z,x).

Using column vectors, matrix multiplication (as defined using linear map composition) is from right to left, so to rotate around fixed axis z (matrix Rz) and then around fixed axis x (matrix Rx), we write the multiplication to get the total rotation matrix R this way: R = Rx * Rz.

Euler angles do not operate on fixed axes, but on moving ones, as shown in the animation.  So we cannot use the fixed axes matrices, at least not directly.
The solution is to reverse the order. This means to perform an Euler angle rotation around z and then (moving axis) x we use: R = Rz * Rx.
But why does this work out? Here is a proof using basic linear algebra.

Proof:
To perform  an Euler angle rotation around z and then (moving axis) x' the matrix composition goes this way: R = Rx' * Rz, where x' is the moving axis, obtained by rotating x by Rz.

We call the basis (or frame) obtained by rotating the standard basis E3 by Rz the basis F'. Using a change of basis, we notice that in the basis F', the rotation around x' is simply a rotation around the 'fixed' axis x (of that basis).

Formally, with rx' being the linear map represented by Rx' in the standard basis, we have this change of basis:

Rep{E3,E3}(rx') = Rep{F',E3}(id) * Rep{F',F'}(rx') * Rep{E3,F'}(id)

Notice that the representation of the basis F' with respect to E3, is Actually Rz, since the column vectors of Rz are the representations in R3 of the rotated E3 vectors, so Rep{F',E3}(id) = Rz and, since rotation matrices are non-singular (invertible) Rep{E3,F'}(id) = (Rz)^-1 .
Additionally, as we explained,  Rep{F',F'}(rx) is simply equal to Rx.

In terms of matrices, this change of basis is: Rx' = Rz * Rx * (Rz)^-1.
Replacing Rx' in the expression of the total rotation, we get:

R
= Rx' * Rz
= ( Rz * Rx * (Rz)^-1 ) * Rz
= Rz * Rx * ( (Rz)^-1 * Rz )
= Rz * Rx * ( I )
= Rz * Rx.

Q.E.D

So, we can always use composition in reverse order (post multiplication) of fixed axis rotation matrices to obtain the total rotation matrix for a sequence of Euler angle rotations.