is used only when M1 and M2 are of comparatively same masses, as in 2 celestial bodies.

well I do not know who said that. the formula F(t) = G * M1 * M2 / r(t)^2

is in fact a universal law of physics not just for gravity force, but also for optics, acoustics, magnetic fields, electric fields and even particle physics.

In the case of a gravitation field as long as the masses can be represented by a point mass, the equation seem valid.

I know there will be many experts from places like stack overflow and gamedevnet that will come back with quantum physics, general gravitation, vending for space and * like that, to those people I can only say this engine does not deal with that.

let us clarify something. in the equation G * M1 * M2 / r(t)^2

we assume that the orbit is mostly circular, so when this is that case the estimation at f at t = t0 + dt

is simply the projection of that position one step in the future to the circulatory trajectory, it works because you have the knowledge that the trajectory is approximated circular. therefore it only change direction but not much in magnitude

for a generic method what to do is that you expand

f (t+dt) = G * M1 * M2 / r(t + td)^2 as an algebraic polynomial expansion and and you discard the higher order term of dt^n

which is what I suspect the expression you are using is.

please try using the term I suggest, see how that goes. it should work as long as your orbits are close to circular. it will also work for elliptic orbits as a long as you know what orbit you want to accomplish.

for any other initial condition value problems, it will not produce the desired result, you need use he expansion but there is not guarantee of a stable orbit.