Behaviour of DC series motors
If I understand you correctly, you have a function \$I\mapsto \omega(I) \$ (the mechanical simulation) and you are looking for \$I\$ that solves
\$ I = { V - G I \omega(I) \over R } \$.

This is a steady state solution, there are no system dynamics involved, so the issues are with the technique you are using to solve the equation.

There is no reason a priori that the 'relaxation' technique will work.

One hack you could try is to update the current \$I\$ incrementally, which 
introduces some numerical 'damping'.

You could try a binary search to find the satisfying value.

You could implement a secant method (essentially estimate the slope by finite differences).