Here is an uncomfortable secret that textbooks tend to whisper: most equations of motion have no formula for the answer. Newton hands you F = ma, you rearrange it into a differential equation for where things go — and then, for almost anything more interesting than a ball thrown in a vacuum, you discover there is no tidy closed-form x(t) to write down. Three planets tugging on each other? No formula. A pendulum swung hard? No elementary formula. Air resistance that grows with speed, a rocket that loses mass as it burns, a charged particle spiralling through a messy magnetic field? No, no, and no.

So how does a spacecraft find Jupiter, or a game make a thousand asteroids drift convincingly, or a climate model run for a century? The trick is to stop asking for the whole future at once and instead take one tiny step at a time. If you know exactly where something is and how fast it's moving right now, and you know the forces on it, then you know which way it's about to head. Nudge it forward by a sliver of time \Delta t, recompute the forces at the new spot, and step again. Repeat a few thousand times and a trajectory unrolls in front of you. That is numerical integration of an ordinary differential equation (ODE): trading one impossible act of algebra for a great many easy acts of arithmetic. This page is about how to take those steps — and, just as importantly, how not to.

Step one: turn F = ma into a first-order system

Newton's law is a second-order equation — it involves an acceleration, which is a second derivative of position:

Every stepping method we'll meet is built to march a first-order equation forward — something of the shape \dfrac{dy}{dt} = f(y, t), "the rate of change of the state is a known function of the state." The universal move that unlocks everything is to promote velocity to a state variable of its own. Instead of one second-order equation for x, keep track of the pair (x, v) and write two first-order equations:

That's it — a second-order law became a first-order system in twice as many variables. In three dimensions the state is six numbers (x, y, z, v_x, v_y, v_z); for N planets it's 6N. The machinery below never needs to know it came from Newton; it just marches the state vector \mathbf{y} forward using its rate of change f(\mathbf{y}, t). Reduce first, then step.

Euler's method: the most honest, most naive step

The simplest possible idea is due to Leonhard Euler in the 1760s and reads almost like common sense: assume the rate of change holds steady across one small step. If you're moving at velocity v_n right now, then over the next \Delta t you'll travel about v_n\,\Delta t; if your acceleration is a_n, your velocity changes by about a_n\,\Delta t. Turn the crank on the whole state:

Where does this come from, precisely? From chopping the Taylor series off after one term. The exact future position is

Euler keeps the straight-line term and discards everything from \Delta t^2 onward. So the error made in a single step — the local truncation error — is of order \Delta t^2. But to integrate over a fixed span of time T you need T/\Delta t steps, and the little errors accumulate. Multiply an error of size \Delta t^2 by a number of steps proportional to 1/\Delta t and you get a global error of order \Delta t. Euler is a first-order method: halve the step size and you (roughly) halve the final error. That is a painfully slow bargain — to get one more decimal digit of accuracy you must take ten times as many steps.

Watch Euler in action — and watch it cheat

Let's integrate the friendliest system in physics: a mass on a spring, the simple harmonic oscillator. The force is F = -kx, so the acceleration is a = -\tfrac{k}{m}x — always pulling back toward the origin. Take k = m = 1, release the mass from x = 1 at rest, and the true motion is a perfect cosine that swings forever between -1 and +1. Its energy E = \tfrac12 m v^2 + \tfrac12 k x^2 should sit rock-steady at 0.5 for all time. Run Euler and read the last column carefully:

The energy climbs, step after step — the oscillation slowly spirals outward instead of tracing a closed loop. Nothing is wrong with the code; this is Euler being Euler. Each straight-line step overshoots the true curved path just slightly outward, and for this system those overshoots always add energy rather than remove it. Left running for a few hundred steps the amplitude balloons absurdly. A method that invents energy from nothing is a serious problem if you were hoping to model, say, a stable planetary orbit.

Seeing the drift: energy over time

The plot below runs both Euler and the fancier RK4 method (coming up next) on that same unit spring and tracks the energy as time marches on. The true energy is the flat line at 0.5. Watch Euler peel away from it while RK4 clings to it.

Same equation, same step size — wildly different faithfulness. The gap is not about doing more steps; it's about being smarter within each step. That is exactly what Runge and Kutta figured out.

Runge–Kutta 4: sample the slope, then average it wisely

Euler's flaw is that it commits to the slope measured at the start of the step and rides it blindly all the way across, even though the true slope is changing the whole time. The fourth-order Runge–Kutta method (RK4), published around 1900, fixes this by taking four trial peeks at the slope inside a single step and blending them into one superior estimate. Writing the state's rate of change as f(\mathbf{y}, t), one RK4 step from \mathbf{y}_n is:

Then step forward using a weighted average of the four, leaning hardest on the two midpoint estimates:

The weights \tfrac16, \tfrac13, \tfrac13, \tfrac16 are not plucked from the air: they are precisely the numbers that make the method's Taylor expansion agree with the true solution all the way out to the \Delta t^4 term. That buys an enormous prize. RK4's local error is O(\Delta t^5) and its global error is O(\Delta t^4). Halve the step and the error drops by a factor of 2^4 = 16, not a measly factor of 2. RK4 costs four force evaluations per step instead of Euler's one, but for any accuracy you actually care about it wins by a landslide — which is why it is the reliable workhorse behind countless simulators.

Read the energy column: it holds at 0.500 to three decimals across the whole run. Four cheap peeks per step, and the spurious drift essentially vanishes.

A real trajectory: a planet in orbit

Now for something with no closed-form shortcut in general — a body falling around a star under gravity. The acceleration points toward the centre and weakens with the square of the distance:

The state is now four numbers (x, y, v_x, v_y). We'll step it with symplectic Euler — a one-line twist on plain Euler that updates the velocity first and then advances the position using that new velocity. That tiny reordering makes the method conserve energy over the long haul (up to a small bounded wobble), which is exactly what you want for an orbit that should stay closed. Set GM = 1 and launch the body from radius 1 with just the right sideways speed for a circle; the energy should hover near E = \tfrac12 v^2 - \tfrac{GM}{r} = -0.5.

The radius stays near 1 and after roughly one period (t \approx 2\pi \approx 6.28) the body arrives back near its launch point: a genuine bound orbit, computed from nothing but the force law and a few hundred additions. Try plain forward Euler on this same setup and the planet spirals slowly out into deep space — a fake, and a catastrophic one for anything meant to run for millions of steps.