the method of lines
A time-dependent PDE has two kinds of derivatives entangled together — derivatives in space and a derivative in time — and discretizing them both at once (as FTCS does) can be confusing and rigid. The method of lines is a clean divide-and-conquer strategy: discretize ONLY the space derivatives, leaving time continuous, which turns the single PDE into a big system of ordinary differential equations in time — one ODE per spatial grid point. Then hand that ODE system to any good ODE solver you like. The name pictures the solution evolving along vertical 'lines' through each grid point as time runs.
Watch it work on the heat equation u_t = k u_xx. Replace only the spatial second derivative by its central difference, leaving the time derivative alone: at grid point i you get du_i/dt = (k/h^2)(u_{i+1} - 2 u_i + u_{i-1}). Do this at every interior grid point and you have a coupled system of ODEs, written compactly as dU/dt = A U, where U is the vector of all the grid values and A is the discrete spatial operator (the stiffness-like matrix). This is now a standard initial-value ODE problem — and you are free to integrate it with whatever time-stepper suits: forward Euler (which reproduces FTCS, with its stability limit), backward Euler or trapezoidal (giving the implicit, unconditionally stable schemes), or a high-order Runge–Kutta. Separating space from time lets you mix and match the best discretization for each.
The method of lines is the conceptual backbone of most modern PDE software: it cleanly separates the spatial discretization (finite differences, finite elements, or spectral) from the time integration, so you can reuse a mature, adaptive, error-controlled ODE library for the hard part of marching in time. It also makes the connection between PDE stability and ODE stiffness vivid — the spatial operator A typically has a huge spread of eigenvalues (very negative for high spatial frequencies on a fine grid), making the ODE system STIFF, which is precisely why explicit time-steppers demand tiny steps (the CFL story) and implicit ones are worth their cost. The main caveat: the approach fits time-evolution (parabolic, hyperbolic) problems naturally, and is less directly suited to purely elliptic, steady-state problems that have no time direction to march along.
Discretize u_t = u_xx on 100 interior points and you get a 100-dimensional ODE system dU/dt = A U whose matrix A has eigenvalues spanning from near 0 down to about -4/h^2 — a stiffness ratio of thousands. Feed it to an implicit, adaptive ODE solver and the time-stepping is handled robustly without you ever choosing dt by hand.
Discretize space, keep time continuous — a PDE becomes a stiff ODE system.
The spatial operator's eigenvalue spread makes the resulting ODE system STIFF, so naively pairing the method of lines with an explicit solver re-imposes the CFL step limit. Stiffness is why implicit time integrators exist; the method of lines makes that dependence explicit.