3 ...  2 ... 1 ... BOOM! 
That's it! Six weeks of hard work.
A 6-DOF fully non-linear rigid body simulator including variable mass written 
in plain C. All equations were derived from scratch. The code is pretty small 
and works on just plain data-structures and has no external dependencies
making future integration into any framework as plain as day.
The rocket (dimmed yellow line) consist of two stages with their centers
indicated by a yellow dot, respectively. The lower stage, the one closer the
bottom, carries double the weight as the upper stage. Hence, the center of 
mass (white dot) is more towards the lower stage, initially. About 90% of the 
lower stage is propellant fuel getting burned at a given rate once the rocket 
lifts off from the ground at full thrust. At the tail end of the lower stage
there is a nozzle indicated with a dimmed white line showing the amount of 
propellant fuel (white: full, black: empty). 
The rocket starts with a given angular velocity at lift-off time. The nozzle 
is held down such that all the thrust is generated in the upward direction 
lifting the rocket off from the ground. During this time the thruster embedded 
within the lower stage produces thrust by burning fuel at a given rate pushing 
it through the nozzle. This consumption of propellant will reduce the mass of 
the lower stage and will also reduces the weight of the rocket and as such 
increase its acceleration even more. And due to the loss in mass of the lower
stage, the center of mass will shift towards the upper stage of the rocket. 
A few seconds after lift-off, but before burnout, the nozzle suddenly tilts
out producing a disaster. In tilt position the nozzle produces a moment/torque 
on the lower stage of the rocket until burnout. Since the thruster is still 
burning at full thrust, the center of mass keeps moves right to the center of 
the upper stage. During this time the center of mass, which also equals the 
center of rotation of the rocket, moves father and father away from the 
nozzle. The generated angular momentum coupled with the change in inertia of 
the lower stage due to propellant burning, and due to the increase in distance 
of the lower stage to the center of mass, will lead the rocket to tune over. 
This becomes only possible this way by carefully considering the non-linear 
gyroscopic term (wxIw) within the angular momentum equation, which is fully 
implemented here, and, in addition, by including the time rate of change of 
the inertia tensor -- which, according to my knowledge, is something new for 
video games, not even speaking about variable mass / center of mass.
The feature of having a rate of change of mass that influences the whole 
system can even be useful for games outside of any realism. For, one can 
consider the rate of change of mass as a parameter. You may specify a constant
rate or make this rate a function all by itself. Since the rate of change of 
mass influences the mass distribution, and as such the motion of the object 
under consideration, one may produce quite interesting motion patterns. That 
is to say, one may build some mass-rate function that will, as an effect, 
shift the center of mass from one point to another and as such influences the
motion of the object. Can be pretty cool, something to experiment with. 
Next-up: 
(a) Runge-Kutta of order 4 for all equations
(b) Wings / Airfoils