function dy = ang_vels(t,y)
    dy = zeros(3,1);    % a column vector
    dy(1) = 1/Ix* (M1-(Iz-Iy)*y(2)*y(3)) ;
    dy(2) = 1/Ix* (M2-(Ix-Iz)*y(1)*y(3)) ;
    dy(3) = 1/Ix* (M3-(Iy-Ix)*y(1)*y(2)) ;