Touch one end of a metal rod to a flame and, a little while later, the far end grows warm. Drop a bead of ink into a glass of still water and watch the dark blot blur outward until the whole glass is faintly grey. In both cases something — heat, dye — spreads: it flows from where there is a lot of it to where there is little, smoothing every bump as it goes. The law behind this is the heat equation (or diffusion equation), and it is a partial differential equation (PDE) — an equation that ties together how a quantity changes in space and how it changes in time at the same moment:

Here u(x,t) is the temperature at position x along the rod at time t, and \alpha is the diffusivity — how eagerly the material shuffles heat around. In words the equation says: a point heats up (or cools) at a rate set by how curved the temperature profile is right there. A sharp peak (strongly curved, \partial^2 u/\partial x^2 < 0) collapses; a valley fills in. Left alone, every wrinkle irons itself flat.

The trouble is that only a handful of PDEs can be solved with pen and paper, and only for the tidiest shapes and boundaries. The rod in your hand, with its lumpy initial temperature and its clamped ends, is not one of them. So we do what all of computational physics does: we put the field on a grid, replace the derivatives with differences, and march the whole thing forward one small step at a time. This page builds that recipe — the finite-difference method — from scratch, runs it live in your browser, and then confronts the one lesson that trips up every newcomer: do it carelessly and the numbers do not just get inaccurate, they explode.

Step 1 — put the field on a grid

A function u(x,t) holds a value at infinitely many points; a computer can only store finitely many numbers. So we sample. Chop the rod into equal slices of width \Delta x and chop time into equal ticks of length \Delta t. We only ever track the temperature at the crossing points of this lattice, and we give each one a pair of indices:

Read it aloud: u_i^{\,n} is "the temperature at the i-th point in space, at the n-th moment in time." The subscript i walks along the rod; the superscript n counts the clock ticks (it is an index, not a power). A whole row u_0^n, u_1^n, \dots, u_M^n is a single snapshot of the rod at one instant. Our job is to turn each snapshot into the next one.

Step 2 — turn derivatives into differences

A derivative is a limit of a ratio of small differences. On a grid we simply stop taking the limit and keep the ratio. That is the whole trick. For the time derivative we look one tick ahead — a forward difference:

For the second derivative in space we compare a point with its two neighbours. If the point sits below the average of its neighbours the profile is cupped upward (a valley); if above, it is a peak. The central second difference measures exactly that:

The pattern (\,\text{left} - 2\times\text{middle} + \text{right}\,) is the discrete curvature: it is zero for a straight line, negative at a hump, positive in a dip. It is worth knowing by heart — it is the workhorse of the whole subject.

Step 3 — the FTCS update rule

Now substitute both approximations into the heat equation \partial u/\partial t = \alpha\,\partial^2 u/\partial x^2 and solve for the one thing we do not yet know — the temperature at the next tick, u_i^{\,n+1}:

Multiply through by \Delta t, gather the constant into a single dimensionless number

and out drops the explicit FTCS update — "Forward-Time, Central-Space":

And that is the entire algorithm: from one snapshot compute the next by sweeping this formula across the interior points, then repeat. We only need two more ingredients to actually run it — where to start, and what to do at the ends.

Starting values and the ends of the rod

A PDE on its own has infinitely many solutions; to pick out our rod we pin down two things.

• Initial condition — the temperature everywhere at t = 0, i.e. the whole starting row u_i^{\,0}. Maybe a hot spike in the middle, maybe a step (left half hot, right half cold).
• Boundary conditions — what happens at the two ends, every tick. The simplest are Dirichlet conditions: hold the end temperatures fixed (say both ends clamped in ice at u = 0). In code we simply never update u_0 and u_M — the sweep runs only over the interior points i = 1,\dots,M-1.

With those in hand the loop is trivial to write. Here it is — a hot spike planted in the centre of an ice-clamped rod, marched forward with a stable r = 0.4. Each printed row is a snapshot: the bars are the temperature of each cell, so you can literally watch the spike slump and spread outward while its peak sinks. Press Run ▶, then change a number and run it again.

The catch: stability

You might reasonably expect that a smaller time step \Delta t always makes the answer better, and that a bigger one just makes it a bit rougher. For this explicit scheme that is false, and the failure is spectacular. Push \Delta t too far and the solution does not degrade gently — it erupts into wild, alternating-sign oscillations that double in size every few steps until the numbers overflow. There is a hard line:

Where does the magic \tfrac12 come from? Here is the intuition (a von Neumann stability argument, stripped to its bones). Any error on the grid can be built from wavy pieces. Feed the single wiggliest wave the grid can hold — one that flips sign at every point, \dots, +1, -1, +1, -1, \dots — into the update. For that pattern u_{i+1} = u_{i-1} = -u_i, so

Each step just multiplies that wave by the amplification factor g = 1 - 4r. The wave stays under control only if |g| \le 1, i.e. -1 \le 1 - 4r \le 1, which rearranges to 0 \le r \le \tfrac12. The moment r > \tfrac12 we get g < -1: every step the jagged wave grows and flips sign — precisely the exploding zig-zag you are about to see. This is a member of the famous CFL family of conditions (Courant–Friedrichs–Lewy) that govern nearly every explicit scheme in computational physics: the time step is chained to the grid spacing, and you cannot cheat it.

Watch it happen. The next block is the same solver as before — same rod, same spike — but with \Delta t nudged up so that r = 0.8 > \tfrac12. We print the raw numbers so you can see the sign flip and the magnitude run away:

Pinning down the threshold

The r = \tfrac12 boundary is sharp — worth seeing directly rather than taking on faith. This last program starts every run from the worst possible initial condition (the sign-flipping zig-zag, the fastest-growing mode) and sweeps a range of time steps. For each it marches 200 steps and reports the peak size. Below the line the peak stays near 1; a hair above it, the numbers detonate:

Everything with r \le 0.5 stays tame; the very first value above it (r = 0.6) has already grown by dozens of orders of magnitude. That is the stability condition, made concrete.

Watching heat diffuse

Run the stable scheme long enough on a concentrated blob of heat and the profile relaxes into a smooth, ever-widening bell that flattens as it spreads — the exact quantity of heat is conserved, it just gets shared out. The graph below shows that idealised spreading: the temperature profile u(x) at four successive times, narrow-and-tall giving way to broad-and-low.

Notice how the tallest, most sharply curved peak changes fastest: strong curvature means a large \partial^2 u/\partial x^2, and the heat equation turns curvature directly into rate of change. As the profile flattens, its curvature drops, and the spreading slows — which is why the last two snapshots look far more alike than the first two.