mjbWorld

This page
http://www.euclideanspace.com/maths/differential/matrixcalculus/
directed me here.
It's hard to tell what exactly the problem is in your first computations, since it's possible that you define things with some different convention than I am used to. But as far as I can tell you are doing things with normal conventions. So I'll try to explain what the right way is as I understand it, and maybe you can get that to work with your system.

Basically, the relation I know of is that Rdot = -omega*R (not Rdot = omega*R). This is for an R, which is the rotation matrix from some frame 1 to frame 2. So a vector in frame 1 would be rotated to frame 2 like this: V2 = R*V1
The rotation matrix is defined by a series of Euler rotation angles. This link shows the way I'm used to defining things:
http://mathworld.wolfram.com/RotationMatrix.html
(equation 5-7) Some of the signs there differ from your Euler page
http://www.euclideanspace.com/maths/geometry/rotations/euler/index.htm
So that may be part of the problem.
For example, where R1 is first computed, the signs on the sine entries are backwards. Your R1 would be the rotation from frame 2 to frame 1 (since it's the transpose of the "normal" way). A handy way of checking your rotation matrices is to say "I want to rotate the frame 1 y-axis by 90 degrees" and then the result should lie along one of the axes of frame 2, which is easy to check with just 1's and 0's in your R1. Check this for each of the axes, and you should be good. If it doesn't do like you expect, something is probably wrong in the rotation matrix.

I'll use x1,y1,z1 for the frame 1 axes and x2,y2,z2 for the frame 2 axes. So looking at the figure above your R1 equation, if you start with the y1 axis and rotate an attitude of 90 deg, the new y2 is along the old -x1 and the new x2 is along the y1. So with that arrangement, a vector in frame 1 that lies along the y1 axis should be the same a vector in frame 2 along the x2 axis (same sign). But if you multiply that out using 90 deg in your R1, you get the following:
[0 -1 0
1 0 0
0 0 1]
So for V2 = R1*V1, if V1 is a unit vector on the y-axis like I was testing above [0 1 0]', then you'll get X2 = [-1 0 0]'. Which is wrong. As I described above a 90 deg rotation of y1 should give you +x1 (not -x1). It's possible that somehow you are doing things with some different sign convention or something, but I don't think you are based on how you are using things. You do seem to use R1 to go from frame 1 to frame 2 like I would think and you use right-hand axes and rotations, so that all seems right.

The other part is omega then, which is the skew-symmetric form of the rotation rates.
So for x,y,z axes, the rate of rotation about those axes is wx, wy, and wz. And the skew-symmetric matrix is
[0 -wz wy
wz 0 -wx
-wy wx 0]
Using those definitions, Rdot = -omega*R should be true.

Hope that helps some and isn't too confusing. If you haven't been there, check out the gamedev.net math and physics forums - they're good about helping with these kinds of things.