Pretty much all of the game development I have been involved with runs afoul of simulating a physical world in discrete time steps. This is of course very simple, but hardly elegant (not to mention mathematically inaccurate). It also has severe disadvantages when large values are involved (either very large speeds, or very large time intervals).

I'm trying to make a continuous physics simulation, just for learning, which goes like this:

time = get_time()
while true do
    new_time = get_time()
    update_world(new_time - time)
    render()
    time = new_time
end

And update_world() is a continuous physical simulation. Meaning that for example, for an accelerated object, instead of doing

object.x = object.x + object.vx * timestep
object.vx = object.vx + object.ax * timestep -- timestep is fixed

I'm doing something like

object.x = object.x + object.vx * deltatime + object.ax * ((deltatime ^ 2) / 2)
object.vx = object.vx + object.ax * deltatime

However, I'm having a hard time with the numerical stability of my solutions, especially for very large time intervals (think of simulating a physical world for hundreds of thousands of virtual years). Depending on the framerate, I get wildly different solutions.

How can I improve the numerical stability of my continuous physical simulations?