Chapter 24: Partial Differential Equations — Finite Differences for Heat, Wave, and Flow

"The purpose of computing is insight, not numbers." — Richard W. Hamming, Numerical Methods for Scientists and Engineers

Overview

For twenty-three chapters your heat solver has been a promise. You chose the problem in Chapter 1, grew a program heat that compiled, made the temperature field a 2D allocatable array in Chapter 5, lifted the update into a step procedure in Chapter 6, split it into modules, and bundled its state into a field_t derived type in Chapter 9. Through all of it, the actual physics inside step has been a deliberate placeholder — a stand-in that moved heat around plausibly enough to lock in the interface while we built everything else. We told you, more than once, that the real thing arrives in Chapter 24.

This is Chapter 24. Here the placeholder becomes a genuine, correct, finite-difference solver for the heat equation — the computational core the entire book has been building toward, and the piece the Chapter 38 capstone will optimize, parallelize, validate, and present as a paper. When you finish this chapter, ./heat will simulate real diffusion on a real grid, and you will understand every line of why it works — including the one condition that, if you get it wrong, makes the whole simulation explode into a blizzard of infinities.

A partial differential equation is how physics writes down a field that changes in both space and time, and finite differences are the oldest, most transparent way to make a computer solve one. The method is almost embarrassingly simple: replace every derivative with a difference of neighbouring grid values, and march. The subtlety — and this is the intellectual heart of the chapter — is that "march" hides a trap. Take a time step even slightly too large and your simulation does not merely lose accuracy; it detonates, doubling its error every step until the numbers overflow. Understanding exactly where that cliff is, and staying on the safe side of it, is what separates a working solver from a random-number generator.

In this chapter, you will learn to:

• Write the heat/diffusion equation in 1D and 2D and the wave equation, and read each term as a physical statement.
• Turn the continuous Laplacian $\nabla^2 u$ into the discrete five-point stencil on a grid, and see why the $1/\Delta x^2$ scaling is not optional.
• March the solution forward with an explicit scheme (FTCS), and understand what an implicit scheme would buy you and what it would cost.
• Find the stability cliff — the CFL condition — derive the 2D limit $r \le 1/4$, compute a safe timestep, and watch a run blow up when you cross it.
• Hold the edges of your grid fixed (Dirichlet), insulated (Neumann), or wrapped (periodic), and code each one correctly.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — this is your chapter. Read every section; §24.2 (the stencil), §24.4 (stability), and the Project Checkpoint are the whole method. Hand-trace the tiny grids yourself — that is where the understanding lives. - 📖 Standard — read straight through. This is the climax of the numerical-methods part and the moment the running project becomes a real simulation. - 🔧 Legacy — the relaxation kernels you modernized in Part IV are exactly the five-point stencil of §24.2, usually with a Gauss–Seidel or Jacobi sweep. Read §24.2 and §24.5 to see what that old code was really computing. - ⚡ HPC — the stencil update in §24.3 and the Project Checkpoint is the kernel you will vectorize in Chapter 29 and parallelize in Chapters 33–34. Read the ⚡ Performance Notes closely; the loop order and the halo you will need later are decided by the physics here.

24.1 From Physics to a Grid: the Heat and Wave Equations

Every simulation in this book — and most of the simulations that run on the world's supercomputers — starts with a partial differential equation. You have already met its cousin. In Chapter 23 an ordinary differential equation described a quantity that changes with respect to one variable, usually time: a planet's position $\mathbf{r}(t)$, a reactant's concentration $c(t)$. But temperature on a metal plate is not one number that changes in time; it is a number at every point of the plate that changes in time. It varies in space and time at once, and to describe how, you need derivatives with respect to more than one variable.

Definition (partial differential equation). A partial differential equation (PDE) is an equation relating a function of several variables to its partial derivatives — its rates of change with respect to each variable separately, written $\frac{\partial u}{\partial t}$, $\frac{\partial u}{\partial x}$, and so on, where the $\partial$ ("partial") signals that the other variables are held fixed. Where an ODE governs a function of one variable, a PDE governs a field: a quantity spread over space that also evolves in time. The heat equation, the wave equation, the equations of fluid flow, electromagnetism, and quantum mechanics are all PDEs.

The heat equation. Picture the square plate from Chapter 1. Let $u(x, y, t)$ be its temperature at position $(x, y)$ and time $t$. Heat flows from hot to cold, and the rate at which a point heats up or cools down is proportional to how much hotter or colder its immediate surroundings are — how curved the temperature profile is around it. That physical statement is the heat equation. In two dimensions:

$$ \frac{\partial u}{\partial t} = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) = \alpha \nabla^2 u $$

The left side is the rate of temperature change at a point. The right side measures the spatial curvature of the temperature field, bundled into the Laplacian operator $\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}$, scaled by the thermal diffusivity $\alpha$ (how readily the material conducts heat — large for copper, small for wood). The one-dimensional version, a heated rod, is just the same idea with one spatial term:

$$ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} $$

💡 Intuition: the Laplacian is a "neighbourhood comparison." At any point it asks: is the average of my neighbours higher or lower than me? If your neighbours are hotter on balance, $\nabla^2 u > 0$ and you warm up; if you are a local hot spot, $\nabla^2 u < 0$ and you cool down; if you are exactly the average of your surroundings, $\nabla^2 u = 0$ and you sit still. That is diffusion in one sentence: everything drifts toward the average of what surrounds it. When nothing is changing anywhere — $\partial u/\partial t = 0$ everywhere — you have reached steady state, and the equation collapses to $\nabla^2 u = 0$, Laplace's equation, the smooth temperature distribution the plate settles into.

The same equation, with $u$ reinterpreted, governs an enormous range of physics: it is the diffusion equation for a chemical concentration spreading through a solvent, for a pollutant dispersing in groundwater, for the probability density of a random walker. "Heat equation" and "diffusion equation" are two names for the identical mathematics, which is why a solver you build for a warm plate transfers, almost unchanged, to a dozen other sciences. This is the first recurring theme in a numerical dress: Fortran fits its domain like a key fits a lock, and its domain is precisely fields-on-grids evolving in time.

The wave equation. Not everything diffuses. Pluck a guitar string, drop a stone in a pond, send a seismic pulse through rock, and the disturbance does not smear out and settle — it propagates, travelling at a definite speed and, in the ideal case, keeping its shape. That behaviour comes from a second time derivative:

$$ \frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u $$

where $c$ is the wave speed. The two-derivatives-in-time structure is what gives waves their momentum: a point that is displaced does not simply relax back, it overshoots, because its velocity carries it through. The heat equation forgets; the wave equation remembers. We will focus the chapter's engineering on the heat equation — it is the project's spine — but the wave equation shares the same spatial Laplacian and the same finite-difference machinery, and its stability condition in §24.4 is the original one that gave the whole subject its name. We will return to it there.

🔗 Connection: in Chapter 23 you learned the method of lines: turn a PDE into a large system of ODEs by discretising space and leaving time continuous, then hand the system to an ODE integrator. That is exactly what this chapter does. Once §24.2 replaces $\nabla^2 u$ with an algebraic expression in the grid values, the heat equation becomes one ODE per grid point — $\frac{du_{ij}}{dt} = \alpha\,(\text{stuff involving neighbours})$ — and §24.3's explicit step is nothing but the forward-Euler method from Chapter 23 applied to that system. You already know the time integrator; this chapter is about discretising the space.

🔄 Check Your Understanding. 1. What makes the heat equation a partial differential equation rather than an ordinary one? 2. In plain words, what does the Laplacian $\nabla^2 u$ measure at a point? 3. What single structural difference between the heat equation and the wave equation makes one diffuse and the other propagate?

Answers

1. The unknown $u(x, y, t)$ depends on more than one variable, so the equation involves partial derivatives with respect to each ($\partial/\partial t$ and $\partial^2/\partial x^2$, …). An ODE's unknown depends on one variable only. 2. It compares a point to the average of its neighbours: positive when the surroundings are hotter on balance (so the point warms), negative at a local hot spot (so it cools), zero when the point equals its neighbourhood average. 3. The order of the time derivative. Heat has a first time derivative ($\partial u/\partial t$) and relaxes toward equilibrium; the wave equation has a second ($\partial^2 u/\partial t^2$), which acts like inertia and makes displacements overshoot and travel.

24.2 Finite-Difference Discretization: the Five-Point Stencil

A computer cannot store a temperature at every point of a continuous plate — there are infinitely many. So we do the only thing we can: lay down a grid, keep the temperature at the grid points only, and approximate the derivatives from those samples. This is the finite-difference method, and its central move you already met, for a single derivative, in Chapter 22. Here we assemble those difference approximations into the full Laplacian.

Lay a uniform grid over the plate: points $(x_i, y_j)$ with spacing $\Delta x$ in the $x$-direction and $\Delta y$ in the $y$-direction. Write $u_{i,j}$ for the temperature at grid point $(i, j)$ — in Fortran, exactly the array element u(i,j). We need the two second derivatives at each interior point.

Recall from Chapter 22 the central second difference, the three-point approximation to a second derivative:

$$ \frac{\partial^2 u}{\partial x^2}\bigg|_{i,j} \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta x^2} $$

Read it as curvature: it is the point's two $x$-neighbours, added, minus twice the point, divided by the spacing squared — larger when the point sits in a valley or on a peak of the profile. This approximation is second-order accurate: its error shrinks like $O(\Delta x^2)$, so halving the grid spacing quarters the error. (Chapter 22 derived that from a Taylor expansion; we will verify it numerically for the whole stencil in the exercises, which is the project increment Chapter 22 asked you to set up.) The $y$-derivative is the identical idea in the other index:

$$ \frac{\partial^2 u}{\partial y^2}\bigg|_{i,j} \approx \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta y^2} $$

Add them, and the continuous Laplacian becomes a formula in five grid values — the point and its four nearest neighbours:

$$ \nabla^2 u \big|_{i,j} \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{\Delta x^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{\Delta y^2} $$

When the grid is square — $\Delta x = \Delta y = h$, the common case — this collapses to something you should commit to memory:

$$ \nabla^2 u \big|_{i,j} \approx \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4\,u_{i,j}}{h^2} $$

Definition (five-point stencil). The five-point stencil is the finite-difference approximation to the 2D Laplacian that combines a grid point with its four nearest neighbours — north, south, east, west — as (sum of the four neighbours) minus (four times the centre), all divided by $h^2$. The five points it touches give it its name. Drawn on the grid, it is a small plus sign:

text (i, j+1) | (i-1, j) -- (i,j) -- (i+1, j) north+south+east+west - 4*centre | ------------------------------- (i, j-1) h^2

It is the discrete workhorse of computational physics: the same plus-sign pattern appears in heat and diffusion codes, in the pressure solves of fluid dynamics, in electrostatics, and in image processing (where it is an edge detector). Learn it once and you recognise it everywhere.

Look closely at that neighbour-sum-minus-four-times-centre — you have seen it before. In the Chapter 5 Project Checkpoint you wrote exactly this expression as a single whole-array statement with array sections, but deliberately without the $1/h^2$. We told you then that the scaling was deferred to Chapter 24. This is Chapter 24, and here is why the scaling is not a detail you can drop.

🚪 Threshold Concept — the $1/h^2$ is the difference between a number and a derivative. The bare sum $u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j}$ is just a comparison of neighbours — a pure number, with no physics in it. Dividing by $h^2$ is what converts that comparison into an actual second derivative, a rate of curvature with real units (temperature per length-squared). Get the scaling wrong — forget it, or use $h$ instead of $h^2$ — and every term downstream is off by a factor that depends on how finely you gridded, so your simulation's answer would change when you refined the mesh, which is the signature of a broken discretisation. The neighbour sum is Chapter 5's array exercise; the $1/h^2$ is the physics. Once you see that a finite-difference stencil is a scaled comparison of neighbours, every PDE discretisation you ever meet is a variation on this one idea.

Here is the stencil computing a Laplacian we can check by hand. Take the smooth field $u(x,y) = x^2 + y^2$, whose exact Laplacian is $\frac{\partial^2}{\partial x^2}(x^2+y^2) + \frac{\partial^2}{\partial y^2}(x^2+y^2) = 2 + 2 = 4$ everywhere. A good discrete Laplacian should return $4$ at every interior point — and because the five-point stencil is exact for quadratics (its error involves fourth derivatives, which vanish here), it returns $4$ exactly, for any spacing:

program laplacian_check
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 5
  real(dp), parameter :: h = 0.5_dp
  real(dp) :: u(n,n), lap(n,n)
  integer  :: i, j

  do j = 1, n
    do i = 1, n                                   ! inner loop over first index (Ch.5 column-major)
      u(i,j) = (real(i-1,dp)*h)**2 + (real(j-1,dp)*h)**2   ! u = x^2 + y^2
    end do
  end do

  lap = 0.0_dp
  lap(2:n-1,2:n-1) = ( u(1:n-2,2:n-1) + u(3:n,2:n-1)      &   ! north + south
                     + u(2:n-1,1:n-2) + u(2:n-1,3:n)      &   ! east  + west
                     - 4.0_dp*u(2:n-1,2:n-1) ) / h**2         ! - 4*centre, scaled

  print '(a)', 'discrete Laplacian of x^2 + y^2 (interior):'
  do i = 2, n-1
    print '(3f8.2)', lap(i,2:n-1)
  end do
end program laplacian_check

$ gfortran -std=f2018 -Wall example-01-five-point-stencil.f90 -o lap && ./lap
discrete Laplacian of x^2 + y^2 (interior):
    4.00    4.00    4.00
    4.00    4.00    4.00
    4.00    4.00    4.00

Verify one point by hand so you trust the machine. At $(2,2)$, with $h = 0.5$: the centre is $u_{2,2} = (0.5)^2 + (0.5)^2 = 0.5$; the $x$-neighbours are $u_{1,2} = 0 + 0.25 = 0.25$ and $u_{3,2} = 1.0 + 0.25 = 1.25$; the $y$-neighbours are $u_{2,1} = 0.25$ and $u_{2,3} = 1.25$. The stencil is $(0.25 + 1.25 + 0.25 + 1.25 - 4\cdot 0.5)/0.25 = (3.0 - 2.0)/0.25 = 1.0/0.25 = 4$. Exactly four, as the mathematics demanded. That single statement — array sections for the neighbours, one division for the scaling — is the spatial heart of your solver.

⚡ Performance Note: written as array sections, u(1:n-2, 2:n-1) + u(3:n, 2:n-1) + …, the stencil hands the compiler four whole-array reads and one whole-array write with no aliasing between input and output — precisely the shape Chapter 27 shows a Fortran compiler vectorising into SIMD instructions that process several grid points per clock. The equivalent nested do loop compiles to the same or faster code as long as the inner loop runs over the first index (Chapter 5's column-major rule), because then consecutive iterations touch consecutive memory. Loop over the second index instead and you stride through memory the long way; on a large grid that single mistake costs you a factor of several, and we will measure it in Chapter 27. The stencil is a memory-bound kernel — it does little arithmetic per value fetched — so how you walk memory is the whole performance game.

🐍 Python Comparison: in NumPy you would write this stencil with array slicing that looks almost identical — lap[1:-1,1:-1] = (u[:-2,1:-1] + u[2:,1:-1] + … - 4*u[1:-1,1:-1]) / h**2 — and for one big vectorised stencil NumPy is genuinely fast, because it is calling down into compiled loops. The gap opens when you must iterate: a heat simulation applies this stencil tens of thousands of times, each step depending on the last, and if any part of that time loop drops into interpreted Python it crawls. The Fortran version keeps the entire time loop compiled. This is the split the book keeps returning to — Python to orchestrate and plot, Fortran for the kernel that runs ten thousand times — and it is why Chapter 15 wrapped exactly this step with f2py.

24.3 Marching in Time: Explicit vs Implicit Stepping

The stencil gives us the spatial side, $\alpha \nabla^2 u$, at one instant. To get a simulation we must also advance time: knowing the temperature field now, compute it a moment $\Delta t$ later, and repeat. As the method-of-lines connection above promised, this is an ODE-integration problem, and the simplest integrator is the forward-Euler method from Chapter 23: approximate the time derivative by a forward difference and evaluate the right-hand side at the current time.

Apply that to the heat equation. The time derivative becomes $\frac{u^{n+1}_{i,j} - u^n_{i,j}}{\Delta t}$, where the superscript $n$ counts time steps ($u^n$ is the field now, $u^{n+1}$ one step later), and the right-hand side is the five-point stencil evaluated on the current field $u^n$:

$$ \frac{u^{n+1}_{i,j} - u^n_{i,j}}{\Delta t} = \alpha \, \frac{u^n_{i+1,j} + u^n_{i-1,j} + u^n_{i,j+1} + u^n_{i,j-1} - 4u^n_{i,j}}{h^2} $$

Every quantity on the right is known — it is the field you already have — so you can solve for the one unknown, $u^{n+1}_{i,j}$, directly. Rearranging gives the update rule that is your solver:

$$ u^{n+1}_{i,j} = u^n_{i,j} + r\left(u^n_{i+1,j} + u^n_{i-1,j} + u^n_{i,j+1} + u^n_{i,j-1} - 4u^n_{i,j}\right), \qquad r \equiv \frac{\alpha\,\Delta t}{h^2} $$

The dimensionless group $r = \alpha \Delta t / h^2$ collects everything — the physics ($\alpha$), the timestep ($\Delta t$), and the grid ($h$) — into a single number that will dominate the rest of the chapter. It is often called the diffusion number (or mesh Fourier number). Hold onto it; §24.4 is entirely about how big it is allowed to be.

Definition (explicit scheme). A time-stepping scheme is explicit when the new value at each point is given by a formula in terms of already-known (current-step) values only — you compute $u^{n+1}$ directly, point by point, with no equations to solve. The scheme above, forward-Euler in time on a centred-space stencil, is the classic explicit method for the heat equation; it is abbreviated FTCS, for Forward-Time, Centred-Space. Explicit schemes are trivial to program (they are just the update formula, swept over the grid) and cheap per step — but, as we are about to see, they buy that simplicity with a hard limit on how large $\Delta t$ may be.

Definition (implicit scheme). A scheme is implicit when the update formula involves the unknown new values $u^{n+1}$ on both sides — for example, evaluating the stencil at the new time level rather than the old. Then $u^{n+1}$ is not given by a formula; it is the solution of a coupled system of linear equations, one per grid point, that you must solve every step. Implicit schemes cost far more per step (a linear solve) but are typically unconditionally stable — they place no limit on $\Delta t$ — so they win decisively when the explicit stability limit would force absurdly tiny steps.

That trade-off — cheap-but-limited versus expensive-but-unconditional — is one of the genuine forks in computational physics, and you are already equipped for both branches. The explicit branch is this chapter. The implicit branch leads straight to the linear algebra of Chapter 21: an implicit heat step is precisely a system $A\mathbf{u}^{n+1} = \mathbf{u}^n$, and you would solve it by calling LAPACK's dgesv, exactly as that chapter's optional Project Checkpoint sketched. For the plate and timescales in this book, explicit stepping is the right tool — simpler, faster to write, and fast enough — so we build it in full and leave the implicit variant as the road not taken (but clearly signposted).

Here is the FTCS update as a self-contained program you can hand-check: a 1D rod, five nodes, ends held at $0$ and $100$ degrees (fixed-temperature boundaries), interior starting cold, marched three steps with a safe $r = 0.25$. The 1D update drops the $y$-neighbours: $u^{n+1}_i = u^n_i + r(u^n_{i+1} - 2u^n_i + u^n_{i-1})$.

program heat_1d
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 5, nsteps = 3
  real(dp), parameter :: r = 0.25_dp            ! r = alpha*dt/dx^2, <= 1/2 in 1D: stable
  real(dp) :: u(n), u_new(n)
  integer  :: i, step

  u      = 0.0_dp
  u(1)   = 0.0_dp                                ! left end held at 0    (Dirichlet)
  u(n)   = 100.0_dp                              ! right end held at 100 (Dirichlet)
  u_new  = u

  do step = 1, nsteps
    do i = 2, n-1                                ! interior only; ends never change
      u_new(i) = u(i) + r*(u(i+1) - 2.0_dp*u(i) + u(i-1))
    end do
    u = u_new                                    ! commit the whole new field at once
    print '(a,i0,a,5f10.4)', 'step ', step, ': ', u
  end do
end program heat_1d

$ gfortran -std=f2018 -Wall example-02-ftcs-1d.f90 -o heat1d && ./heat1d
step 1:     0.0000    0.0000    0.0000   25.0000  100.0000
step 2:     0.0000    0.0000    6.2500   37.5000  100.0000
step 3:     0.0000    1.5625   12.5000   45.3125  100.0000

Trace step 1 to convince yourself. Only interior nodes $i=2,3,4$ update. Nodes $2$ and $3$ see all-zero neighbourhoods, so they stay $0$. Node $4$ sees the hot end: $u_{\text{new}}(4) = 0 + 0.25\,(100 - 2\cdot 0 + 0) = 0.25 \times 100 = 25$. Step 2, node 4: $25 + 0.25\,(100 - 50 + 0) = 25 + 12.5 = 37.5$; node 3, now with a warm neighbour: $0 + 0.25\,(25 - 0 + 0) = 6.25$. The heat is visibly crawling leftward from the hot end, one node per step, and every value the program prints is a dyadic fraction you can reproduce with a pencil. Left to run, this rod settles to the straight line $0, 25, 50, 75, 100$ — the 1D steady state, where the temperature profile is linear and $\nabla^2 u = 0$ everywhere.

⚠️ Common Pitfall — update from a snapshot, never in place. Notice the code writes into a separate u_new and only copies it back after the whole sweep. This is not fussiness; it is correctness. The FTCS formula says every new value depends on the old neighbours, $u^n$. If you wrote back into u in the same loop, then by the time you computed node $i$, node $i-1$ would already hold its new value, and you would be mixing time levels — computing a different scheme (in fact a Gauss–Seidel relaxation) without meaning to. Sometimes that is what you want; here it is a silent bug that quietly changes your physics. Two buffers, swapped each step, keep the time levels clean. (The whole-array step in the checkpoint gets this for free: it reads the field into the Laplacian before touching it.)

🔄 Check Your Understanding. 1. In the FTCS update $u^{n+1}_{i,j} = u^n_{i,j} + r(\dots)$, why can you solve for $u^{n+1}$ directly, with no linear system — and what makes that "explicit"? 2. What is the diffusion number $r$, in terms of $\alpha$, $\Delta t$, and $h$? 3. Why does the code update into u_new and copy back, rather than overwrite u as it sweeps?

Answers

1. Every term on the right-hand side is evaluated at the current time level $n$, which is fully known, so the new value is a plain arithmetic expression in known quantities — that "known-values-only" property is exactly what explicit means. An implicit scheme would put unknown $u^{n+1}$ values on the right too, coupling all points into a system you must solve. 2. $r = \alpha\,\Delta t / h^2$ — physics times timestep over spacing squared; the dimensionless group that controls both accuracy and stability. 3. Because FTCS is defined in terms of the old neighbours $u^n$. Overwriting in place would feed already-updated neighbours into later points, silently switching to a different (Gauss–Seidel) scheme.

24.4 Stability: the CFL Condition

Now the trap. You have a correct stencil and a correct update rule, and you might reasonably think that a smaller timestep is merely slower and a larger one merely less accurate — that $\Delta t$ trades speed for precision and nothing worse. For the heat equation solved explicitly, that intuition is dangerously wrong. There is a hard threshold in $\Delta t$. Below it the simulation is stable and marches sensibly toward steady state. Above it — even a hair above it — the simulation does not just lose accuracy; it becomes violently unstable, with errors that double, or worse, every single step until the numbers overflow to infinity. This threshold is the single most important thing to understand in this chapter, and one of the most important in all of computational physics.

Where does the cliff come from? Watch what the FTCS update does to a single interior point surrounded by its own opposite — the worst case for the scheme. Concretely, consider the "checkerboard" pattern, where adjacent grid points alternate high and low, $u_{i,j} = A\,(-1)^{i+j}$. This is the most rapidly varying field the grid can represent, and diffusion should crush it fast. Feed it through the update: each of the four neighbours has the opposite sign to the centre, so the neighbour sum is $-4$ times the centre, and

$$ u^{n+1}_{i,j} = u^n_{i,j} + r\,(-4u^n_{i,j} - 4u^n_{i,j}) = (1 - 8r)\,u^n_{i,j}. $$

Every step multiplies this pattern's amplitude by the factor $G = 1 - 8r$. For the simulation to be stable, that factor must not grow the pattern: we need $|G| \le 1$, i.e. $|1 - 8r| \le 1$. The upper side is automatic; the binding constraint is the lower side, $1 - 8r \ge -1$, which rearranges to $8r \le 2$, or

$$ \boxed{\,r = \frac{\alpha\,\Delta t}{h^2} \le \frac{1}{4}\,} \qquad\Longleftrightarrow\qquad \Delta t \le \frac{h^2}{4\alpha}\quad(\text{2D, } \Delta x=\Delta y=h). $$

Cross that line — take $r > 1/4$ — and $G = 1 - 8r < -1$, so the checkerboard pattern flips sign and grows every step. Any real field contains a whisper of this pattern (floating-point round-off alone seeds it), and once $|G| > 1$ that whisper is amplified geometrically until it swamps everything and overflows. That is the detonation. A fuller von Neumann stability analysis — feeding a general wave $e^{\mathrm i(k_x x + k_y y)}$ through the scheme instead of just the checkerboard — gives the amplification factor $G = 1 - 4r[\sin^2(k_x h/2) + \sin^2(k_y h/2)]$, and the checkerboard ($k_x h = k_y h = \pi$, both sines $=1$) is exactly the worst mode, confirming $r \le 1/4$.

Definition (CFL condition; stability). A numerical scheme is stable when errors present in the solution (from round-off, from the initial data) stay bounded as it marches, rather than growing without limit. For an explicit scheme, stability requires the timestep to be small enough — a restriction named the CFL condition, after Courant, Friedrichs, and Lewy, who identified it in 1928. The physical idea: in one timestep, the numerical scheme only "sees" as far as its stencil reaches (one grid cell), so $\Delta t$ must be small enough that whatever the physics does in that time also stays within about one cell. For the heat equation the concrete form is the diffusion-number limit $r = \alpha\Delta t/h^2 \le 1/4$ in 2D (and $\le 1/2$ in 1D, $\le 1/6$ in 3D — the pattern is $1/(2d)$ in $d$ dimensions).

⚠️ A precise word on the name. Strictly, "CFL condition" was born for wave (hyperbolic) problems, where information travels at a finite speed $c$ and the condition is the Courant number $C = c\,\Delta t/h \le 1$ — literally "a wave may not cross more than one cell per step." The heat equation is parabolic; its explicit-stability limit $r \le 1/4$ is a diffusion-number (von Neumann) condition, not a Courant number. In everyday practice, though, computational scientists say "the CFL condition" loosely for any explicit scheme's timestep-stability limit, and this book follows that common usage — while being honest that the two are different beasts. The practical distinction matters enormously: the wave limit scales like $\Delta t \sim h$, but the diffusion limit scales like $\Delta t \sim h^2$. Halve your grid spacing to resolve finer detail, and an explicit heat solver demands a four-times-smaller timestep (and thus four times as many steps) — the notorious reason explicit diffusion codes get painfully slow on fine grids, and the reason implicit methods exist.

Let us make the cliff visible. The program below is the same FTCS solver, but run unstably: a $4\times 4$ plate with fixed (zero) edges, seeded with a checkerboard in its four interior cells, and stepped with $r = 0.5$ — double the stable limit. Watch the interior.

program cfl_blowup
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 4
  real(dp), parameter :: r = 0.5_dp             ! r = 0.5 > 1/4 : UNSTABLE in 2D
  real(dp) :: u(n,n), u_new(n,n)
  integer  :: i, j, step

  u = 0.0_dp                                     ! zero Dirichlet edges, held fixed
  u(2,2) =  1.0_dp;  u(2,3) = -1.0_dp            ! checkerboard on the 2x2 interior
  u(3,2) = -1.0_dp;  u(3,3) =  1.0_dp

  call show(0, u)
  do step = 1, 3
    u_new = u
    do j = 2, n-1
      do i = 2, n-1
        u_new(i,j) = u(i,j) + r*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1) - 4.0_dp*u(i,j))
      end do
    end do
    u = u_new
    call show(step, u)
  end do

contains
  subroutine show(s, a)
    integer,  intent(in) :: s
    real(dp), intent(in) :: a(:,:)
    integer :: ii
    print '(a,i0,a)', 'step ', s, ':'
    do ii = 1, size(a,1);  print '(4f9.2)', a(ii,:);  end do
  end subroutine show
end program cfl_blowup

$ gfortran -std=f2018 -Wall example-03-cfl-blowup.f90 -o blow && ./blow
step 0:
     0.00     0.00     0.00     0.00
     0.00     1.00    -1.00     0.00
     0.00    -1.00     1.00     0.00
     0.00     0.00     0.00     0.00
step 1:
     0.00     0.00     0.00     0.00
     0.00    -2.00     2.00     0.00
     0.00     2.00    -2.00     0.00
     0.00     0.00     0.00     0.00
step 2:
     0.00     0.00     0.00     0.00
     0.00     4.00    -4.00     0.00
     0.00    -4.00     4.00     0.00
     0.00     0.00     0.00     0.00
step 3:
     0.00     0.00     0.00     0.00
     0.00    -8.00     8.00     0.00
     0.00     8.00    -8.00     0.00
     0.00     0.00     0.00     0.00

Hand-trace the first step at $(2,2)$: it starts at $1$, its two interior neighbours $(2,3)$ and $(3,2)$ are each $-1$, and its two boundary neighbours are $0$, so $u_{\text{new}}(2,2) = 1 + 0.5\,(0 + (-1) + 0 + (-1) - 4\cdot 1) = 1 + 0.5\,(-6) = 1 - 3 = -2$. The amplitude went from $1$ to $2$ and the sign flipped. Every step repeats the trick: $1 \to -2 \to 4 \to -8$, a factor of $-2$ each time, magnitude doubling without bound. Continue the run and the interior reaches $\pm 16, \pm 32, \dots$; after about a thousand steps the values exceed $\sim 1.8\times 10^{308}$, the largest real(dp) — the number Chapter 20 called huge(1.0_dp) — and overflow to Infinity, and one more step of Inf - Inf gives NaN. The simulation is destroyed.

🐛 Find the Bug. A colleague's heat code "was working yesterday" but now fills the screen with NaN after a few hundred steps. The only change: they refined the grid from $101\times 101$ to $201\times 201$ points "for a sharper picture," halving $h$, and left dt and alpha untouched. What happened, and what is the one-line fix?

Answer

Halving $h$ quadruples $r = \alpha\Delta t/h^2$, because $r$ depends on $h^2$. Their old run must have had $r$ between $1/16$ and $1/4$; the finer grid pushed $4r$ past $1/4$, across the CFL cliff, and the scheme went unstable. The fix is to recompute the timestep from the grid every time the grid changes: dt = 0.9_dp * h**2 / (4.0_dp*alpha) (a 10% safety margin under the limit), never hard-code it. This is the single most common way real explicit diffusion codes blow up, and why the project's step will always derive dt from the grid, not trust a constant.

The lesson is a rule you never break with an explicit scheme: compute the timestep from the stability limit, do not guess it. The Project Checkpoint below builds that rule into a small function, and the capstone's validation study leans on it. Choose $r$ safely under $1/4$ — a common choice is $r \approx 0.2$ or a fixed fraction like $0.9 \times \frac14$ — and the solver is rock-solid; drift over $1/4$ and no amount of clever coding saves you.

And the wave equation? Its explicit scheme (central differences in both space and time) carries the original CFL condition in its pure form: $C = c\,\Delta t / h \le 1$ in 1D. Here the interpretation is crisp and physical — the numerical wave must not outrun the true wave, so a disturbance may travel at most one grid cell per timestep. It is the same principle as the heat limit (an explicit stencil sees only one cell per step), wearing the clothes it was originally cut for.

📜 From History: the condition is named for a 1928 paper by Richard Courant, Kurt Friedrichs, and Hans Lewy — "Über die partiellen Differenzengleichungen der mathematischen Physik" ("On the partial difference equations of mathematical physics"). Remarkably, they were not trying to run simulations; there were no electronic computers. They were using finite differences as a theoretical tool to prove that certain PDEs have solutions at all, and they discovered that the difference scheme only converges to the true solution if the time and space steps are suitably related. Two decades later, when von Neumann and others actually programmed these schemes on the first machines, that abstract convergence condition turned out to be the concrete line between a simulation that works and one that explodes. It has governed explicit PDE codes ever since — a rule discovered on paper a generation before the machines it would rule.

24.5 Boundary Conditions: Dirichlet, Neumann, Periodic

The stencil at an interior point reaches to its four neighbours. But what about a point on the very edge of the grid, which is missing neighbours off the boundary? A PDE alone does not determine its own solution; you must also say what happens at the edges of the domain. These are the boundary conditions, and they are as much a part of the physics as the equation itself — the same heat equation on the same plate gives completely different answers depending on whether the edges are held hot, insulated, or wrapped around. Three kinds cover the overwhelming majority of problems.

Definition (Dirichlet boundary condition). A Dirichlet condition fixes the value of the field on the boundary: $u = g$ on the edge, for some prescribed $g$. For the plate, "hold the top edge at 100°C and the others at 0°C" is Dirichlet. In code it is the simplest of all: set the boundary grid points to their prescribed values once (or re-impose them each step if they vary in time), and never update them — the time loop sweeps only the interior. Every example so far in this chapter used Dirichlet edges.

Definition (Neumann boundary condition). A Neumann condition fixes the normal derivative of the field on the boundary: $\partial u/\partial n = q$. The most common case is $q = 0$, an insulated (or "adiabatic," "no-flux," "zero-gradient") edge: no heat crosses it. You impose it discretely by forcing the field to have zero slope across the edge — the boundary point is set equal to its inward neighbour, u(1,j) = u(2,j), so the difference across the edge, and hence the flux, is zero. (A "ghost cell" just outside the domain, mirrored from just inside, is the tidy general way to do this.)

Definition (periodic boundary condition). A periodic condition wraps the domain so that the last point's neighbour is the first point: the field is treated as tiling space seamlessly, $u_{n+1} \equiv u_1$. Use it to model a small representative patch of a large uniform system without edge effects — a ring, a crystal lattice, a slice of atmosphere. In code the neighbour index is taken modulo the grid size, so the stencil at the last point reaches around to the first. Fortran's modulo intrinsic (Chapter 3) gives the correct non-negative wrap.

The three are worth seeing side by side as code. Suppose u(1) is a left edge in a 1D field of length n; here is how each condition sets it (Dirichlet and Neumann shown; periodic couples the two ends):

! Dirichlet: the edge holds a fixed value, set once and never updated in the sweep
u(1) = 0.0_dp                       ! left end pinned at 0 degrees

! Neumann (zero-flux / insulated): edge mirrors its inward neighbour -> zero gradient
u(1) = u(2)                         ! no heat crosses the left edge

! Periodic: the domain wraps; index arithmetic is done modulo n on the interior update
!   left neighbour of point i is  1 + modulo(i-2, n),  right is  1 + modulo(i, n)
im = 1 + modulo(i-2, n)             ! wraps 1 -> n
ip = 1 + modulo(i,   n)             ! wraps n -> 1
u_new(i) = u(i) + r*(u(ip) - 2.0_dp*u(i) + u(im))

💡 Intuition: the three conditions answer three different physical questions about the edge. Dirichlet: "what is the temperature here?" — a wall clamped to a thermostat. Neumann: "how much heat flows across here?" — most often "none," a perfectly insulated wall. Periodic: "what is on the other side?" — the other side of the same domain, as if the plate were rolled into a cylinder. Choosing wrongly does not crash the code; it silently solves a different problem. A steady-state plate with one hot Dirichlet edge and three cold ones reaches a smooth temperature gradient; make all four edges insulated (Neumann) instead and, with no heat entering or leaving, the whole plate drifts to a single uniform temperature — the average of where it started.

For the running project we use Dirichlet everywhere: fixed-temperature edges on the square plate, which matches the Chapter 1 problem statement ("hot on one edge, cold on the others") and keeps the solver's core clean. Neumann and periodic edges are the natural first extensions — the exercises and case studies build them — and the domain-decomposition of Chapter 34 turns inter-process boundaries into a fourth kind, the "halo" or ghost-cell exchange, which is periodic coupling between one process's edge and its neighbour's.

⚠️ Common Pitfall — the corners, and the off-by-one at the edge. Two boundary mistakes recur. First, loop bounds: an $n \times n$ grid with a one-cell boundary updates only i = 2, n-1 and j = 2, n-1. Write do i = 1, n by reflex and the stencil reads u(0,j) or u(n+1,j) — out of bounds, an instant crash with -fcheck=all (Chapter 13) or, without it, silent memory corruption. Second, corners: a corner cell belongs to two edges, so if the top edge is $100$ and the left edge is $0$, the top-left corner is over-determined. It rarely matters for Dirichlet (corners do not enter any interior stencil), but decide deliberately rather than letting statement order decide for you.

🔄 Check Your Understanding. 1. You want to model a metal bar whose two ends are clamped to ice baths at $0°C$. Which boundary condition, and what does the update loop do to the end points? 2. What is the physical meaning of a zero-Neumann edge, and how is it coded in one line? 3. Why does a periodic boundary use modulo rather than a plain mod or a hand-written if?

Answers

1. Dirichlet. Set both end values to $0$ once and exclude them from the sweep (do i = 2, n-1); they are held fixed, never updated. 2. Zero flux — a perfectly insulated edge that no heat crosses. Code it by setting the edge equal to its inward neighbour, e.g. u(1) = u(2), forcing zero gradient across the boundary. 3. modulo(i-2, n) returns a result with the sign of the divisor (always non-negative here), so the wrap at both ends is correct; mod can return a negative index at the low end, and a hand-written if is easy to get subtly wrong. modulo is the one-intrinsic answer.

24.6 Structured Grids and Output for Visualization

Everything in this chapter has lived on a structured grid, and it is worth naming the choice, because it is the reason the code is so clean.

Definition (structured grid). A structured grid is a mesh whose points are arranged in a regular array, so that a point's neighbours are found by simple index arithmetic — $(i\pm1, j)$ and $(i, j\pm1)$ — with no connectivity information stored. Your u(i,j) is the structured grid: the array's own layout encodes who neighbours whom. The alternative, an unstructured grid (arbitrary triangles or polyhedra, as used for complex geometries like an aircraft), must store an explicit list of each cell's neighbours and is far more involved. Structured grids map perfectly onto Fortran arrays, which is a large part of why Fortran and finite-difference methods grew up together.

The regularity is a gift that keeps giving. It is why the Laplacian is a handful of array sections rather than a neighbour-lookup; it is why the update vectorises (Chapter 27); and it is why the plate decomposes so naturally across processors — you just cut the array into rectangular blocks (Chapters 32–34). A great deal of high-performance scientific computing is structured-grid computing precisely because the structure buys all of this.

Finally, a simulation you cannot see is hard to trust. Your solver produces a 2D array of temperatures at each step; to watch heat spread you write those arrays to disk and view them. The full treatment — industry-standard VTK files that open in ParaView, complete with a time-series animation of your plate warming — is Chapter 26, and the project increment there wires it in. For now, the simplest possible output is a plain text dump of the field, which gnuplot or a three-line Python script with matplotlib.pcolormesh will render as a heat map. Even a text print of a small grid, as in every example above, is a form of visualisation — you can see the hot edge and the warm cells creeping inward. Keep the output routine (write_field, from your Chapter 7 I/O work) decoupled from the solver, and you can upgrade it from text to VTK later without touching the physics.

🔗 Connection: decoupling output from computation is not just tidy; it is how real codes stay portable across a decade of changing file formats. A climate model's dynamical core and its I/O layer are separate modules for exactly this reason — the science outlives the storage format. Your heat_io module (Chapter 8) is a small rehearsal of that discipline, and Chapter 25 will give it a self-describing format (NetCDF/HDF5) for runs too large for text.

Project Checkpoint

This is the one you have waited twenty-three chapters for. We now replace the placeholder physics inside step with the real finite-difference heat solver: the five-point stencil, an explicit FTCS update, a CFL-safe timestep, and Dirichlet boundaries — all operating on the canonical field_t from Chapter 9. The public signature step(field, alpha, dt), frozen since Chapter 6, does not change; only its body becomes real, and — as Chapter 9 promised — field is now a type(field_t), so the routine finally has the grid spacing field%dx it needs for the $1/h^2$ scaling.

Two procedures go into heat_solver, matching the canonical interfaces (_style-bible.md §4): a pure function laplacian(u, dx, dy) that returns the properly scaled Laplacian, and the step that uses it. A small helper, stable_dt, computes a CFL-safe timestep from the grid so no one ever hard-codes it:

pure function laplacian(u, dx, dy) result(lap)
  real(dp), intent(in) :: u(:,:), dx, dy
  real(dp) :: lap(size(u,1), size(u,2))
  integer  :: nx, ny
  nx = size(u,1);  ny = size(u,2)
  lap = 0.0_dp                                        ! boundaries: Laplacian unused (held fixed)
  lap(2:nx-1,2:ny-1) =                                              &
        ( u(1:nx-2,2:ny-1) - 2.0_dp*u(2:nx-1,2:ny-1) + u(3:nx,2:ny-1) ) / dx**2   &
      + ( u(2:nx-1,1:ny-2) - 2.0_dp*u(2:nx-1,2:ny-1) + u(2:nx-1,3:ny) ) / dy**2
end function laplacian

subroutine step(field, alpha, dt)                     ! frozen signature; field is now field_t
  type(field_t), intent(inout) :: field
  real(dp),      intent(in)    :: alpha, dt
  real(dp), allocatable :: lap(:,:)
  integer :: nx, ny
  nx = field%nx;  ny = field%ny
  lap = laplacian(field%u, field%dx, field%dy)        ! evaluate on the OLD field
  field%u(2:nx-1,2:ny-1) = field%u(2:nx-1,2:ny-1)  &  ! forward-Euler update, interior only
                         + alpha*dt * lap(2:nx-1,2:ny-1)   ! Dirichlet edges untouched
end subroutine step

pure function stable_dt(alpha, dx, dy, safety) result(dt)
  real(dp), intent(in)           :: alpha, dx, dy
  real(dp), intent(in), optional :: safety            ! fraction of the CFL limit (default 0.9)
  real(dp) :: dt, s
  s = 0.9_dp;  if (present(safety)) s = safety
  dt = s / (2.0_dp*alpha*(1.0_dp/dx**2 + 1.0_dp/dy**2))   ! 2D limit; = s*h^2/(4 alpha) if dx=dy
end function stable_dt

Three design choices earn their place. The Laplacian is a pure function — no side effects — which both documents intent and lets the Chapter 27/29 optimiser reason freely about it. The step evaluates laplacian on the current field before writing anything, so it reads $u^n$ cleanly and the snapshot correctness of §24.3 is automatic (no old = field copy needed). And stable_dt encodes the §24.4 limit once, using the general $\Delta x \neq \Delta y$ form $\Delta t \le 1/[2\alpha(1/\Delta x^2 + 1/\Delta y^2)]$, which reduces to the memorable $h^2/(4\alpha)$ on a square grid.

The full self-contained program — bundling small kinds, heat_types, and heat_solver modules with a driver — is code/project-checkpoint.f90. It builds a $5\times 5$ plate with $\Delta x = \Delta y = 1$, holds the top edge at $100°$, sets $\alpha = 1$ and a safe $\Delta t = 0.2$ (so $r = 0.2 \le \frac14$), and takes two steps. Hand-computed, the interior after step 1 is a single warm row, and after step 2 the warmth has spread inward and sideways, symmetric about the centre column:

after step 1:                          after step 2:
  100.00 100.00 100.00 100.00 100.00     100.00 100.00 100.00 100.00 100.00
    0.00  20.00  20.00  20.00   0.00       0.00  28.00  32.00  28.00   0.00
    0.00   0.00   0.00   0.00   0.00       0.00   4.00   4.00   4.00   0.00
    0.00   0.00   0.00   0.00   0.00       0.00   0.00   0.00   0.00   0.00
    0.00   0.00   0.00   0.00   0.00       0.00   0.00   0.00   0.00   0.00

Trace the centre cell $(2,3)$ across both steps. At step 1 its neighbours are still the initial field — $100$ above and three $0$s — so $0 + 0.2\,(100 + 0 + 0 + 0 - 4\cdot 0) = 0 + 0.2(100) = 20$. At step 2 it now sees $100$ above, $0$ below, and two warm $20$s beside it: $20 + 0.2\,(100 + 0 + 20 + 20 - 4\cdot 20) = 20 + 0.2(60) = 32$ — the warmest interior cell, as the left–right symmetry demands. Compare this to the placeholder output from Chapters 6 and 8 (0, 10, 10, 0 in the second row): that stand-in used $\alpha\Delta t = 0.1$ and no $1/h^2$ — it was, in hindsight, this very stencil with $h = 1$ and no stability check. The solver was always shaped like the real thing; this chapter gives it the correct scaling, a principled timestep, and the physics to back it.

This is the core the rest of the book optimises and scales. Chapter 26 makes it write VTK so you can watch it; Chapter 27 measures its loop order; Chapter 29 tunes the stencil; Chapters 33–34 parallelise the update behind this same step interface; and the Chapter 38 capstone validates it against the analytical solution and presents it as a paper. From here on, your solver solves.

Summary

This chapter turned the running project from a plausible placeholder into a correct, stable finite-difference solver for the heat equation — the computational heart of the book.

The two things to memorize. First, the five-point stencil with its $1/h^2$: the neighbour sum minus four times the centre, divided by the spacing squared — that is the discrete Laplacian, and it is exact for quadratics. Second, and above all, the stability limit $r = \alpha\Delta t/h^2 \le 1/4$ for the explicit 2D scheme: it is the line between a simulation and an explosion, it scales as $\Delta t \sim h^2$, and you compute your timestep from it — you never guess.

Spaced Review

Retrieval practice on the two chapters this one is built from: arrays (Chapter 5) and ODEs (Chapter 23). Answer before peeking.

1. (Ch. 5) The stencil update sweeps u(i,j) in a nested loop. Recalling column-major order, which index should the inner loop run over for cache-friendly access, and why?
   

   Answer

   The **first** index, `i`. Fortran stores arrays column-major, so `u(i,j)` and `u(i+1,j)` are adjacent in memory; running the inner loop over `i` walks memory contiguously. Loop over `j` inside instead and every iteration jumps a whole column — the classic slowdown Chapter 27 measures on exactly this kernel.

2. (Ch. 5) Why can the whole interior Laplacian be written as one array-section statement, u(1:n-2,2:n-1) + u(3:n,2:n-1) + …, instead of a double loop — and what does the compiler gain?
   

   Answer

   Because array **sections** are first-class arrays: each shifted slice is itself an array, and whole-array `+` acts elementwise. The one statement hands the compiler four non-overlapping array reads and one write with no aliasing, which it can vectorise — *more* information than a loop, which is why Fortran's array syntax is fast, not just short.

3. (Ch. 23) The explicit heat step is forward-Euler applied to a system of ODEs. What is that reframing called, and what plays the role of the ODE right-hand side $f$?
   

   Answer

   The **method of lines**: discretise space, leave time continuous, and you get one ODE per grid point, $du_{ij}/dt = f$, where $f = \alpha\nabla^2 u$ — the (scaled) five-point stencil evaluated at that point. FTCS is forward-Euler on that system.

4. (Ch. 23) Euler's method has a local error of $O(\Delta t^2)$ per step. Given that, and the stencil's $O(h^2)$ spatial error, what limits the accuracy of the FTCS heat solver, and how would Chapter 23's toolkit improve the time side?
   

   Answer

   Both contribute: time error $O(\Delta t)$ globally (forward-Euler is first-order in time) and space error $O(h^2)$. The time accuracy is the weaker link; swapping forward-Euler for a higher-order integrator (RK4, or Crank–Nicolson for diffusion) raises the time order — the method-of-lines view makes that a drop-in change of integrator, exactly Chapter 23's point.

5. (Ch. 23) Chapter 23 warned that stiff ODEs force implicit methods. In what sense is the discretised heat equation stiff, and how does that connect to §24.4?
   

   Answer

   The system's fastest-decaying modes (the checkerboard) relax on a timescale $\sim h^2/\alpha$, far shorter than the physical timescale of interest — the hallmark of stiffness. An explicit method must resolve that fastest mode for *stability* (the $r \le 1/4$ limit), forcing tiny steps; an implicit (unconditionally stable) method escapes the restriction — the classic stiff-equation trade-off, here wearing a PDE's clothes.

What's Next

Your solver now computes real physics, and it can print a small grid you can read at a glance. But a $1000 \times 1000$ plate stepped ten thousand times is ten billion numbers no one will ever read as text, and a text file that large is both enormous and slow. Chapter 25 opens Part VI with the self-describing binary formats — NetCDF and HDF5 — that real scientific codes use to store fields at scale, portably and with their metadata attached. Then Chapter 26 makes your solver write VTK, so you can finally open ParaView and watch the heat spread across your plate, one timestep at a time — the payoff image of the whole project. The physics is done; now we make it visible.