Hi Ross,
You have solved the first four equations specified by Earles correctly as far as I can see. However, although he states they are coupled equations, they aren't! Not really. They are coupled only in the sense that they are all subject to a common forcing function(F), but there is no feedback from one to another. For example, q1 doesn't depend in any way on q2, q2dot, q3, q3dot,...etc. You could solve each pair (q and qdot) independently and get the same results. It is only when Earles makes his approximations, and, in effect, adds extra relationships that there is any interaction and potential for instability.
Have you tried solving his second set of equations? You will probably find some instabilities there.
Edit:
On looking at this more carefully, I note that F is not an external forcing function, but merely indicates the transfer of forces between the disc and the other parts. I think you have to use the approximation that q4 = q1 - r.sin(theta).q3 etc. and eliminate q4 and F from the equations, leaving you with six equations (3 for dqi/dt and 3 for dqidot/dt). You need to rearrange these as suggested in my first reply so that you don't have more than one dqidot/dt in any one equation. Then you must introduce a perturbation by, for example, specifying an initial displacement for q1(0) say. Depending on the size of the perturbation you might or might not see an instability.
Edit 2
I've had a go at doing the above myself. Needless to say, what I've done needs very careful checking, but it does seem to show a boundary between stable and unstable regimes.
Alan
