A puzzle: the right answer, the wrong-looking step
Here is a small scene that has confused generations of newcomers. You are solving the initial-value problem y'(t) = -1000*(y - cos(t)) - sin(t), starting from y(0) = 0, and the true solution is gentle: after a brief blink it tracks cos(t) almost exactly, a smooth wave that never changes faster than cos does. So you reach for the trusty workhorse of guide 1, forward Euler, pick a step like h = 0.01 that would resolve cos(t) beautifully, and press go. The result is not a small error — it is a catastrophe. The numbers swing to a hundred, then ten thousand, then overflow to infinity within a dozen steps. The solution you are chasing is calm, yet the method explodes.
Nothing is wrong with your code, and nothing is wrong with the solution. The blow-up is not the accuracy failure you met in guides 2 and 3, where a too-large h gave a recognisable, bounded error you could shrink. This is a different beast: an instability that has nothing to do with how smooth the answer is. The equation contains a hidden fast scale — that factor of -1000 — that decayed away in the first instant and is now physically irrelevant, yet it still dictates the step size your explicit method is allowed to take. An equation that does this has a name. It is stiff, and stiffness is the single most important reason this subject keeps an entire family of implicit methods in its toolbox.
The test equation and the region of absolute stability
To see exactly what goes wrong, we strip the problem to its bones. Linearise any decaying behaviour and you get the test equation y' = lambda*y, where lambda is (often) a negative real number — the true solution is y(t) = y_0 * exp(lambda*t), which decays to zero. The number lambda*h, step size times rate, is what every method really cares about. Apply forward Euler to the test equation: y_{n+1} = y_n + h*lambda*y_n = (1 + h*lambda)*y_n. So each step multiplies the value by the factor (1 + h*lambda). The exact solution shrinks every step, so for the numerical solution to behave we need that factor to have magnitude at most 1: |1 + h*lambda| <= 1.
The set of values of z = h*lambda (allowing z to be complex, since real systems have oscillating modes too) for which a method's growth factor stays inside the unit disk is the method's region of absolute stability. For forward Euler the condition |1 + z| <= 1 is a disk of radius 1 centred at -1: a small, bounded region clinging to the left of the origin. Now plug in our culprit. With lambda = -1000 and h = 0.01, we have h*lambda = -10, and |1 + (-10)| = 9, far outside the disk. Each step multiplies the error by 9 — that is the explosion, exactly. To get back inside we need h*lambda >= -2, i.e. h <= 2/1000 = 0.002, forced by a transient that died long ago.
Going implicit: backward Euler and the price it pays
Here is the rescue. Instead of evaluating the slope at the point you are leaving, evaluate it at the point you are arriving at — that is the backward Euler method: y_{n+1} = y_n + h*f(t_{n+1}, y_{n+1}). Notice the unknown y_{n+1} appears on both sides; the formula is implicit, which is the whole explicit-versus-implicit distinction in one line. Apply it to the test equation y' = lambda*y: y_{n+1} = y_n + h*lambda*y_{n+1}, which rearranges to y_{n+1} = y_n / (1 - h*lambda). The growth factor is now 1/(1 - h*lambda).
Look at what this factor does for our stiff case. With h*lambda = -10, the factor is 1/(1 - (-10)) = 1/11, comfortably less than 1 — the step decays, no matter how big the step. In fact, for every lambda with a negative real part, |1/(1 - h*lambda)| <= 1 for all h > 0. The region of absolute stability is the entire left half of the complex plane. A method whose stability region contains the whole left half-plane is called A-stable, and A-stability is the gold standard for stiff problems: it means the stability limit on the step vanishes entirely, and you may choose h purely for accuracy. Our exploding problem now runs happily at h = 0.01, or even h = 0.1.
Nothing is free. Implicitness has a price: at every step you must solve for y_{n+1} rather than just compute it. For our linear test equation that was easy algebra, but for a real nonlinear system f(t, y) it means solving a (usually nonlinear) algebraic system each step — and the tool for that is exactly Newton's method from the root-finding rung, generalised to systems. So a stiff solver hides a Newton iteration inside every time step, and each Newton step solves a linear system involving the Jacobian of f. That is real work; an implicit step can cost many times an explicit one. The trade is almost always worth it: you replace thousands of forbidden tiny explicit steps with a handful of expensive-but-allowed large implicit ones.
What stiffness really is — and how to recognise it
Now we can give stiffness an honest definition. A stiff equation is one in which the solution contains decay (or oscillation) rates spread over wildly different scales — a fast transient that vanishes almost instantly alongside a slow component you actually want to follow. For a system y' = f(t, y), the relevant rates are the eigenvalues of the Jacobian matrix; stiffness shows up as a large ratio between the most-negative and least-negative real parts, the stiffness ratio. Our scalar example faked a system: the -1000 was the fast mode, cos(t) the slow one, a ratio of about a thousand. Chemical kinetics, electrical circuits with both fast and slow components, and the spatial discretisation of a diffusion PDE are classic stiff sources.
A common misconception is worth puncturing: stiffness is not the same as the solution changing fast. After the initial transient our solution is as smooth and slow as can be, yet the problem is still stiff. The fast scale is present in the equation even when it is absent from the solution — it governs how a perturbation would decay, and an explicit method must respect it at every step regardless. That is the cruel irony of stiffness: explicit methods are forced to take steps tiny enough to resolve a mode that the true solution has long since forgotten. A useful working symptom: if forward Euler or RK4 demands an absurdly small h to stay bounded, far smaller than the smoothness of the answer would suggest, you are almost certainly looking at a stiff problem.
BDF, higher order, and an honest caveat
Backward Euler is A-stable but only first order, so it is too inaccurate to be the everyday workhorse for stiff problems. The production answer comes from the multistep family of guide 4: the backward differentiation formulas (BDF). A BDF fits a polynomial through the new point and several past points, then forces its derivative to match f at the new time — implicit, like backward Euler, and built for stiffness. BDF1 is exactly backward Euler; BDF2 reaches second order while staying A-stable; orders up to 6 are used in practice. These are the engines inside famous stiff solvers, and they are why a library routine that auto-detects stiffness will quietly switch over to a BDF when it sees the symptoms.
Step back and see the whole rung in one frame. Guide 1 gave you Euler and the idea of following the slope; guide 2 separated local from global error and pinned down order; guide 3 built RK4 for high accuracy per step; guide 4 added adaptive step-size control and the efficient multistep methods. This guide added the missing axis — stability — and with it the reason the toolbox is plural. Explicit methods (RK4, Adams-Bashforth) are cheap per step and perfect for non-stiff problems; implicit methods (backward Euler, BDF) cost a Newton solve per step but conquer stiffness. The mature practitioner does not pick a favourite; they diagnose the problem, and let the diagnosis pick the method. And it is worth remembering the bedrock honesty of the whole subject: every number a solver returns is an approximation in finite-precision arithmetic, so stability buys you a bounded, trustworthy answer — never an exact one.