First define each of the axis-aligned rotation matrices as a function of some angle:

rx[θ_]:=(({{1, 0, 0}, {0, Cos[θ], -Sin[θ]}, {0, Sin[θ], Cos[θ]}}));
ry[θ_]:=(({{Cos[θ], 0, Sin[θ]}, {0, 1, 0}, {-Sin[θ], 0, Cos[θ]}}));
rz[θ_]:=(({{Cos[θ], -Sin[θ], 0}, {Sin[θ], Cos[θ], 0}, {0, 0, 1}}));

Now, given Euler angles and assuming the rotation order X, Y, Z with fixed axes, we can convert from Euler angles to a matrix with:

eulertomatrix[θx_, θy_, θz_] := 
      rz[θz] . ry[θy] . rx[θx]

We can now ask Mathematica to write out the matrix for an arbitrary set of input angles:

MatrixForm[eulertomatrix[θx, θy, θz]]

Using the notation Cx to mean Cos[θx], Sz for Sin[θz] and so forth, we can makethe above parameterized matrix look a little cleaner:

MatrixForm[
      eulertomatrix[θx, θy, θz] /. {Cos[θx] -> Cx, Cos[θy] -> Cy, Cos[θz] -> Cz,  
          Sin[θx] -> Sx, Sin[θy] -> Sy, 
          Sin[θz] -> Sz}]

Also, show what the matrix looks like for θy= ±90°

Note that this also shows how gimbal lock kills one degree of freedom.

MatrixForm[
      FullSimplify[
        eulertomatrix[θx, θy, θz] /. {Cos[θy] -> 0, Sin[θy] -> 1}]]
MatrixForm[
      FullSimplify[
        eulertomatrix[θx, θy, θz] /. {Cos[θy] -> 0, Sin[θy] -> -1}]]

Now we can use the above to work out formulae for the Euler angles if we have the matrix:

Note that I think ArcTan takes arguments in reverse order from the C function atan2.

matrixtoeuler[m_] := 
      Block[{x, y, z, cy, ϵ = 0.000001}, 
y = 
          ArcSin[-m[[3, 
                1]]]; 
If[(Abs[
                  m[[3, 1]]] - 
                1) < -ϵ, 
z = 
            ArcTan[m[[1, 1]], 
              m[[2, 
                1]]]; 
x = 
            ArcTan[m[[3, 3]], 
              m[[3, 
                2]]];
, 
x = 
            ArcTan[-m[[1, 3]]/
                m[[3, 
                  1]], -m[[1, 
                    2]]/
                m[[3, 
                  1]]]; 
z = 
            0;
]; 
{x, y, 
          z}
]

Test for some random value

tm = eulertomatrix[ .4,  .9, -1.2];
MatrixForm[tm]
{tx, ty, tz} = matrixtoeuler[tm]
MatrixForm[tmm = eulertomatrix[tx, ty, tz]]
Norm[tm - tmm, ∞]

Test for some random value where middle axis is at 90°

tm = 
        eulertomatrix[1.4, π/2, - .9];
MatrixForm[tm]
{tx, ty, tz} = matrixtoeuler[tm]
MatrixForm[tmm = eulertomatrix[tx, ty, tz]]
Norm[tm - tmm, ∞]

tm = 
        eulertomatrix[- .6, -π/
            2, -1.9];
MatrixForm[tm]
{tx, ty, tz} = matrixtoeuler[tm]
MatrixForm[tmm = eulertomatrix[tx, ty, tz]]
Norm[tm - tmm, ∞]

Test for some more values using a set of loops

\(\(Do[
Do[
Do[

tm\  = \ 
              eulertomatrix[θx, θy, θz]; 
{tx, ty, tz} = matrixtoeuler[tm]; 
tmm = 
              eulertomatrix[tx, ty, 
                tz]; 

\ 
            err\  = \ 
              Norm[tm - 
                  tmm, ∞]; 

If[
              err\  > \ 0.0001, \ 
              Print["Bad: (e=", 
                err, ") ", θx, " ", θy, \ " ", \ θz]];

, {θz, 0. , 
              2  π, 2  π/32}]
, {θy, 0. , 
            2  π, 2  π/32}]
, {θx, 0. , 
          2  π, 2  π/32}];\)\)