I think I've figured out the reason. I cheked the codes and found the value of Izz was filled in my solid body's moment of inertia. So yes, that field should be filled in the moment of inertia relative to its CG. But!! This field seems only set for Izz, so if I want to assign a value to Ixx or Iyy (my axis of rotation), I need to make the z-axis of the global coordinate system lie along the axis of rotation. In other words, I initially set a value for Ixx, but the programme gave this value to Izz, 'cause I didn't know it worked in that way. I think this was the main problem that caused the calculation diverged at early iterations and strange phenomenon.