The numbers are rock steady, except for both apogee and perigee, which still drift, but only by about 10 meters, every 3 seconds...

ah that's a long the line of the error I was expecting.

from here you can reduce that erro even further by doing this secund trick.

the equation of motionof a particle under the action of gravity and not otherhe forces is

m * a = f

a = dv/dt = f / m

we know f / m = -G

where G ar any give position and time is give by -G0/(r * (x + h)

there are two ways to integrate that equation: forward euler and backward euler

forward euler is what the engine does for the velocity and backward for the position.

this makes sence because G close to earth is constantt and does not chanage direction no matter how fast you move near earth.

foward euler is expresed by this

(v1 - v0) / dt = G(t,x)

v1 = v0 + G(t,x) * dt

this is G is the derivative at time and location t, and you mutiply by the dt and add that to the current velocity. graphically this is like this.

- Untitled.png (2.71 KiB) Viewed 4428 times

as you can see when g(t,x) is highly non linear or v0 is very high the change of G(t,x) * dt is quite large and as you can see in the picture v1 drift fromr the ideal v1 that is project on the curve.

let us see the backwork euler now.

this is given by using the velocity not at time v0 but at time v(-dt) and project forward by makinge a change of variable, the final equation is this.

v1 = v0 + G(t + dt, x + dx) * dt

this is, using the value of the derivative at one dt in the future as if it was the derivative at time t.

this is far more complex for function of multiple variables because the way G(t +dt, x, dx) is calcuatione is by expanding on a taylor series and negleting the high order terms.

but for the of an or it busy, is quite eassy to estime G at t + dt.

besially we know that not matter how fast a rocket can travel, r does not changes significantly in magnitude so the linear derivative can be ignored.

the direction however does change significanlly at the speed of a satelity, for example an stationary one has 24 houd period, and rocket move even faster.

graphically the backword euler looks like this.

- Untitled1.png (4.46 KiB) Viewed 4428 times

as you can see in backward euler the velocity at time t + dt end up below the ideal curve, so backward Euler is stable because in loses energy. but that picture is exaggerated, in reality is quite acurate and stable for smooth function like your case.

what remind is how to get G(t +dt, x+dx)

but that is quite eassy, your vehicle is in orbist and has a radius, you know that orbit move on plane, therefore you can stimate G(t+dt, x+dx) by taking the position of the satelite a t0, that is the position of the com of the rigid body.

them get the linar velocity of the rigid body, and add the value v * dt

this is let x0 be the position vector

x1 = x0 + v0 * dt

now with that position, you know the radious is r + h

so you calculate a unit vector for position x1 and multiply by the mangity of r + h, that will be the radius of the rocket at g (t+dt, x+dx)

now you use this value of G to calcuate the force that will be applied to the satellite.

do it right and the error should reduce to a very small value. almost zero.