Attached are two matlab files for plotting and crudely animating rocket evolving 6DOF attitude and position values. This uses quaternions, the quaternion evolution equation, euler's rigid body equations, and Newton's second law. Pcon in the files is the hamiltonian conjugate momenta of the quaternions, but they are not actively used.
The number of letters in the variables for the vectors describe how many components are in them, Quat has 4, Ang has 3 (angular velocity), Anm has 3 (angular momentum), Angv has 4 (angular velocity with 0 as first component), Pcon has 4.
Inertia matrix should be symmetric, otherwise, bad things may happen. (it's always symmetric in real life for reasons)