I have this playground where I try to copy the physical body movement from (in this case) a gizmo mesh. The way I test it is go into debug view, select gizmoMesh, enable position or rotation mode for it. Then when I move/rotate the gizmoMesh the physical body follows.
Its a basic proportional control system. Position is easy and stable. I apply negative difference between gizmo and body position as linear velocity and its sorted. There is a problem with rotation.
For rotation I find the angle between the gizmo and the body and apply a negative euler angle to angular velocity. The problem is that for some angles of one axis, other axis angles might flip 360deg which will throw the control out of balance completely. I suppose its maybe related to gimbal lock.
Below is the playground to test it. The screenshot shows one set of angles that makes it unstable (z around 90)

Quick answer: don’t use Euler angles!!!
There are an infinite number of Euler values to represent a single quaternion orientation.
Moreover, interpolation, gimbal lock, order of axis… happen, so, try to use quaternion instead.
Euler should only be present for user interface.

thanks @Cedric. I was thinking there must be some gimbal going on. The thing is that for proportional control system I need to counter act on the velocity (angular in this case). I am not sure if you can do it without moving to euler.
I read this thread about teleport method for the physics body but it just overwrites all position/quaternion and zeros velocities. It sucks because it overwrites all the forces so physics is lost.
The goal is to control the path of the body while keeping the physics if something interacts with body. It works nice with one axis. you can throw it off balance and it comes back to its path.

You can keep track of the physics by getting/storing/applyin linear and angular velocities.
You can do a an orientation difference with Quaternion just like you would with Euler.

I use quaternions in the example. First conjugate and then multiply quaternions to get the difference in orientation.
Angular velocity is a vector3 so I need to apply velocities as 3 x,y,z values. So I thought I have to move to euler to apply angular velocity?

I checked the code and I would have done the same way as you did.
Are you sure angular velocity is in world space?
Maybe it needs to be changed (each frame) to local space.

That was my thought too but from debugging it doesnt look this way. I tried to move to local by rotating the mesh and body by the same amount so I always reference to no rotation and then calc the difference but didnt help.

I tried with Matrix and rotation seems stable. The only thing its the opposite way I thought difference in matrices is transposed then multiply but must be wrong.