The explicit scheme for the heat equation is a delight to program — one line of arithmetic per grid point — but it comes with a leash. As the stability analysis shows, it is only conditionally stable: the mesh ratio must obey r = \alpha\,\Delta t/\Delta x^2 \le \tfrac12. Refine the spatial grid to chase some fine detail — halve \Delta x — and the largest legal time step drops by a factor of four, because the bound scales with \Delta x^2. On a genuinely fine grid you can be forced into millions of absurdly tiny steps just to keep the arithmetic from exploding.

The cure is to stop computing the new values one at a time and instead solve for the whole new time level at once. Schemes that do this are implicit, and the most celebrated member of the family — the workhorse of practical diffusion solvers since the 1940s — is the Crank–Nicolson method. It buys something remarkable: it is unconditionally stable (no cap on \Delta t at all) and, at the same time, second-order accurate in both space and time. The price is a small linear system solved once per step — and, as we will see, that system is tridiagonal, so the cost is almost as cheap as the explicit rule anyway.

One dial: the θ-method

The explicit and implicit schemes are not really two separate ideas — they are the two ends of a single dial. Start from the heat equation u_t = \alpha\,u_{xx} and, as always, replace the time derivative by a forward difference between level n and level n+1. The only question is where to evaluate the spatial second difference

— at the old level, the new level, or somewhere in between? The θ-method refuses to choose and takes a weighted average, with a single knob \theta \in [0,1]:

Turning the dial recovers three classic schemes:

• \theta = 0 — the whole spatial term sits at the old level: the familiar explicit (forward-time) scheme, conditionally stable.
• \theta = 1 — the whole term sits at the new level: the fully implicit (backward-time) scheme, unconditionally stable but only first-order accurate in time.
• \theta = \tfrac12 — an even split between old and new: Crank–Nicolson. Averaging the two one-sided differences is exactly the trapezoidal rule in time, and that symmetry is what lifts the time accuracy to second order.

So Crank–Nicolson is not an exotic new construction — it is simply the θ-method with the dial set dead-centre, averaging the old and new spatial curvatures.

The update: a tridiagonal system

Set \theta = \tfrac12 and write the mesh ratio r = \alpha\,\Delta t/\Delta x^2. Multiplying through by \Delta t and gathering every new-level value on the left, every old-level value on the right, gives the Crank–Nicolson update at each interior point i:

The right-hand side is entirely known — just old numbers combined with arithmetic, exactly like the explicit scheme. The left-hand side is where the difference lives: three unknown new values are tied together in one equation. Writing this equation at every interior grid point stacks up into a linear system

where \mathbf{d} collects the known right-hand sides (plus the boundary values) and A is tridiagonal: the constant 1+r runs down the main diagonal, and -\tfrac{r}{2} sits on the two neighbouring diagonals — nothing anywhere else. That special shape matters enormously: a tridiagonal system of M unknowns is solved by the Thomas algorithm in O(M) work — one forward sweep and one back-substitution — so an implicit step costs only a small constant multiple of an explicit one, not the O(M^3) of a general matrix solve.

Why it is unconditionally stable

Run the same von Neumann test as before: feed in a single Fourier mode u_i^n = G^n e^{\mathrm{i} k x_i} and ask by how much one step multiplies it. Substituting into the Crank–Nicolson update and writing s = \sin^2\!\big(\tfrac{k\,\Delta x}{2}\big), the exponentials cancel and leave the amplification factor

Here is the whole argument in one line. Because r > 0 and s \ge 0, the quantity 2rs is never negative, so the denominator 1 + 2rs is always at least as large in size as the numerator 1 - 2rs:

No condition on r appears anywhere — the inequality holds for r = 0.5, for r = 50, for r = 5000. Contrast the explicit scheme, whose factor G = 1 - 4r s sails past -1 the moment r > \tfrac12. Crank–Nicolson simply cannot blow up. (Compare the fully implicit scheme, G = 1/(1+4rs), which is also unconditionally stable but loses the second-order time accuracy.)

One step, in code

The entire method is: build the known right-hand side from the old row, then solve the tridiagonal system for the new row. Because the matrix is tridiagonal, the solve is the classic Thomas algorithm — a forward sweep to eliminate the sub-diagonal, then a back-substitution. Here is one Crank–Nicolson step on the interior points, with the ends held at fixed (Dirichlet) boundary values:

Notice there is no stability check anywhere — no test that r \le \tfrac12, because none is needed. You are free to pick \Delta t purely for accuracy, which is exactly the freedom the explicit scheme denies you.

Watch the amplification factor

This is the whole stability story on one picture. The horizontal axis is k\,\Delta x, sweeping from the smoothest mode (left) to the jagged highest-frequency mode k\,\Delta x = \pi (right); the vertical axis is the size |G| of the amplification factor. Anything above the dashed |G| = 1 line grows each step — an instability. Drag the mesh ratio r upward: the explicit curve |1 - 4r\,s| lifts its right end clean over the line the instant r passes 0.5, while the Crank–Nicolson curve \big|(1-2r s)/(1+2r s)\big| stays pinned at or below 1 no matter how far you push — visibly, permanently safe.