Don't discretize everything at once
Everything so far in this rung treated space and time on the same footing: you laid a grid in both, you replaced every derivative — the time derivative and the space derivative alike — by a stencil, and out came a single recipe like FTCS or Crank-Nicolson that told you how to get the whole next time level from the current one. That works, but it bundles two very different jobs into one formula, and it forces you to re-derive a fresh scheme for every new time-stepping idea. The method of lines (MOL) is a change of attitude that untangles the two jobs. Its slogan is simple: discretize space, but leave time alone — for now.
Take the heat equation u_t = alpha * u_xx on a rod. Lay down only the spatial grid: points x_1, x_2, ..., x_M spaced h apart, with the boundaries held fixed. At each interior point define u_i(t), the temperature at that fixed location as a function of time. Now do the half-discretization. Leave the time derivative u_t exactly as it is — a genuine derivative in t. But replace the space derivative u_xx by the central second difference stencil, (u_{i-1} - 2*u_i + u_{i+1}) / h^2. The PDE at point x_i becomes du_i/dt = (alpha/h^2) * (u_{i-1} - 2*u_i + u_{i+1}). That is not a PDE any more. It is an ordinary differential equation in t, one for each grid point.
Picture why the name fits. Each interior grid point sits at a fixed x and traces its temperature forward through time — a horizontal line drawn across the space-time plane. The PDE has been sliced into a bundle of these lines, one per node, and along each line the rule is now an ODE. The lines are not independent: the equation at x_i reaches to its neighbors x_{i-1} and x_{i+1}, so the whole bundle is one big coupled system of ODEs that must be advanced together. But that is a problem from the previous rung, not a new one.
A PDE becomes a giant system of ODEs
Stack the unknowns into one vector u(t) = (u_1(t), u_2(t), ..., u_M(t)) — the whole temperature profile of the rod at time t, frozen into a list of numbers. The central-difference stencil acts the same way at every interior point, so the bundle of ODEs collapses into a single tidy matrix equation: du/dt = A u, where A is the familiar tridiagonal second-difference matrix with -2 down the diagonal and +1 on each side, scaled by alpha/h^2 (the discrete Laplacian in one dimension). This is exactly an initial-value problem for a system of ODEs — the same shape as y' = f(t, y) from the ODE rung, just with y a long vector and f the linear map A. The initial condition is the starting temperature profile u(0); the boundary values feed into the first and last equations as known constants.
Here is the payoff, and it is large. You no longer need to invent a time-stepping scheme for the PDE — you reach into the ODE rung and pull out any solver you like. Apply forward Euler to du/dt = A u and you recover exactly the FTCS scheme; that is not a coincidence, it is FTCS's true identity. Apply backward Euler and you get BTCS. Apply the trapezoidal rule — the average of the slopes at the start and end of the step — and you get Crank-Nicolson. The three schemes you met as separate inventions are just three different ODE integrators bolted onto the same spatial discretization. And nothing stops you from bolting on a better one: a Runge-Kutta method, an adaptive solver, a multistep method, all of which you already understand.
PDE: u_t = alpha * u_xx on 0 < x < 1
step 1 (space): grid x_i = i*h, replace u_xx by central stencil
du_i/dt = (alpha/h^2) * (u_{i-1} - 2 u_i + u_{i+1})
stack into a vector u = (u_1, ..., u_M):
du/dt = A u A = (alpha/h^2) * tridiag(1, -2, 1)
step 2 (time): hand du/dt = A u to ANY ODE solver
forward Euler -> FTCS (explicit, cheap, CFL-limited)
backward Euler -> BTCS (implicit, solve A x = b each step)
trapezoidal -> Crank-Nicolson(implicit, 2nd order in time)Why the CFL limit reappears — and why it is really stiffness
MOL does more than save labor: it explains the stability results of the earlier guides from a single vantage point. A linear system du/dt = A u is stable under a given ODE solver exactly when, for every eigenvalue lambda of A, the scaled value h_t * lambda (time step times eigenvalue) lands inside that solver's region of absolute stability — the very region you studied in the ODE rung. Here A is real and symmetric, a symmetric eigenproblem, so all its eigenvalues are real, which keeps the picture on the real axis. So the question 'is my PDE scheme stable?' becomes 'where do the eigenvalues of A fall, and does my time step shrink them into the stability region?' One picture, used twice.
Now do the arithmetic for the heat equation. The eigenvalues of the second-difference matrix A range from near 0 down to about -4*alpha/h^2 — all real and negative, with the most negative one growing like 1/h^2 as you refine the grid. Forward Euler's stability region only reaches a distance 2 along the negative real axis, so stability demands h_t * (4*alpha/h^2) <= 2, that is h_t <= h^2 / (2*alpha). There is the explicit time-step cap of the FTCS guide, falling straight out of the eigenvalues — and the punishing h^2 is no surprise once you see it comes from that 1/h^2 in A. Backward Euler and Crank-Nicolson have stability regions covering the whole left half-plane, so those negative eigenvalues never escape no matter how large h_t is: unconditionally stable, exactly as the implicit-scheme guide claimed.
The honest fine print
MOL is powerful but not magic, and the same disciplines from the rest of this ladder apply. First, the spatial discretization carries its own truncation error: the central second difference is only O(h^2) accurate, so even with a flawless time integrator your answer is no better in space than that stencil allows. Pairing a high-order ODE solver in time with a low-order stencil in space is wasted effort — the accuracy is set by the weaker of the two. Second, every step still runs in floating-point, where 0.1 has no exact binary form, so shrinking h forever does not keep helping: past a point round-off in the 1/h^2 stencil dominates, the same round-off floor that limits any finite-difference derivative.
Third, choosing an implicit time integrator does not make the work vanish — it relocates it. Backward Euler and Crank-Nicolson require solving a linear system A x = b at every time step, and you rarely form a matrix inverse to do it. In one dimension A is tridiagonal, so the Thomas algorithm solves it in O(M) flops, a bargain. In two or three dimensions the discrete Laplacian is sparse but no longer tridiagonal, and you reach for an iterative solver with a good preconditioner, or a multigrid method — the tools from the iterative-methods rung — precisely because forming or inverting the matrix is out of the question at scale.
One more honest caveat. The clean eigenvalue story holds tidily for linear, constant-coefficient problems on simple grids; for nonlinear PDEs or wild coefficients, A becomes a state-dependent map and the analysis is local and approximate, a guide rather than a guarantee. And for purely advective problems like u_t + a*u_x = 0, the semi-discrete eigenvalues are imaginary rather than negative real, which changes which ODE solvers stay stable and reintroduces numerical dispersion — a reminder that MOL organizes the analysis but does not exempt you from doing it.
What you carry out of this rung
Step back and see how the pieces lock together. Earlier guides taught you to lay a grid and turn derivatives into stencils; to weigh explicit against implicit; to test a scheme for consistency, stability, and convergence and trust the Lax equivalence theorem that consistency plus stability equals convergence; and to find the stable step with von Neumann analysis and the CFL condition. The method of lines is the unifying idea that makes all of that feel like one subject rather than four tricks: discretize space to get a system of ODEs, then bring the entire ODE toolbox — solvers, order, stability regions, stiffness — to bear on time.
It also points outward. Finite differences are the gentlest way to discretize space, but not the only one: the next rung's finite-element (Galerkin) and finite volume methods produce a different A — and a mass matrix M on the left, giving M*du/dt = A u — yet the second half of the recipe is unchanged, because once you have a system of ODEs in time you already know what to do. Spectral methods do the same with a denser, more accurate A. The method of lines is the hinge between 'how do I represent a function on a grid?' and 'how do I march it through time?', and almost every serious PDE code in the world is built around that hinge.