First some definitions.  I'll be rotating about the x axis by theta1, about
the y axis by theta2, and about the z axis by theta3.

Let c = cos(theta1)    s = sin(theta1)
    d = cos(theta2)    t = sin(theta2)
    e = cos(theta3)    u = sin(theta3)
    
Let I be the identity matrix
Let x^t denote the transpose of matrix x
Let x^-1 denote the inverse of matrix x

For part one, the general rotation matrix is simply the product of the
three matrices required to rotate about specific axes.  These matrices
are:

    |1  0 0 0|      |d 0 -t 0|      |e  t 0 0|
x = |0  c s 0|  y = |0 1  0 0|  z = |-t e 0 0|
    |0 -s c 0|      |t 0  d 0|      |0  0 1 0|
    |0  0 0 1|      |0 0  0 1|      |0  0 0 1|

                                            |d*e         d*t        -t  0|
The general rotation matrix is then = xyz = |s*t*e+c*-t  s*t*t+c*e  s*d 0|
                                            |c*t*e+-s*-t c*t*t+-s*e c*d 0|
                                            |0           0          0   1|

Now on to part two.

Note that s^2 + c^2 = t^2 + d^2 = u^2 + e^2 = 1

With those trigonometric identities in mind, note that:

       |1 0        0         0|
xx^t = |0 c*c+s*s  c*-s+s*c  0| = I
       |0 -s*c+c*s -s*-s+c*c 0|
       |0 0        0         1|

       |d*d+-t*-t 0 d*t+-t*d 0|
yy^t = |0         1 0        0| = I
       |t*d+d*-t  0 t*t+d*d  0|
       |0         0 0        1|

       |e*e+t*t  e*-t+t*e  0 0|
zz^t = |-t*e+e*t -t*-t+e*e 0 0| = I
       |0        0         1 0|
       |0        0         0 1|

With the above three matrix equations in mind...
(xyz)^t = (xyz)^-1 <==>
      (xyz)(xyz)^t = I <==>             (by def of matrix inverse)
      (xyz)((x^t)(y^t)(z^t)) = I <==>   (since (ab)^t = (b^t*a^t) and associativity of matrix multiplication)
      x(y(zz^t)y^t)x^t = I <==>         (by associativity of matrix multiplication)
      x(yIy^t)x^t = I <==>              (substitution...)
      xIx^t = I <==>
      I = I

Since these are all bidirectional implications, (xyz)^t does indeed
equal (xyz)^-1