In numerical simulation, a system is called stiff when the state contains multiple time scales, and the fastest mode forces a very small numerical step size for stability, even if the slower modes are the ones we actually care about. This is a numerical property of the problem + method combination, not only a property of the physical model itself.

A common sign of stiffness is the following mismatch: the physical state evolves slowly, but an explicit solver requires extremely small steps to avoid instability. In control engineering, this often appears in models with fast parasitic dynamics, highly damped modes, chemical kinetics, thermal-electrical multi-scale couplings, or strongly dissipative mechanical subsystems.

For a linearized system near an operating point, \( \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{b} \), stiffness is often associated with eigenvalues of \( \mathbf{A} \) having very different magnitudes of real parts.

flowchart TD
  A["Dynamic model xdot = f(t, x)"] --> B["Simulate with explicit method"]
  B --> C["Need tiny h for stability?"]
  C -->|yes| D["Likely stiffness (fast modes constrain h)"]
  C -->|no| E["Nonstiff or mild stiffness"]
  D --> F["Use implicit / stiffness-aware solver"]
  F --> G["Solve algebraic step equation"]
  G --> H["Stable larger time steps for slow behavior"]
        

The standard tool for stability analysis is the scalar test equation \( \dot{y} = \lambda y \), where \( \lambda \in \mathbb{C} \). If \( \Re(\lambda) < 0 \), the exact solution decays:

\[ y(t) = y(0)e^{\lambda t}, \qquad |y(t)| = |y(0)|e^{\Re(\lambda)t}. \]

The exact solution is stable (decaying), but a numerical method may still become unstable if the step size \( h \) is too large. This is the key distinction between continuous-time stability and numerical stability.

For linear systems \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), each mode behaves like the scalar test equation along eigen-directions. If \( \lambda_i \) are eigenvalues of \( \mathbf{A} \), then the relevant numerical parameters are \( z_i = h\lambda_i \). The method must be stable for all \( z_i \) reached by the chosen step size.

A practical stiffness indicator (high-level) for stable linear systems is the stiffness ratio:

\[ \kappa_s = \frac{\max_i |\Re(\lambda_i)|}{\min_{i:\Re(\lambda_i)\neq 0} |\Re(\lambda_i)|}. \]

Large \( \kappa_s \) suggests strong time-scale separation. This is not a complete definition (especially for nonlinear/time-varying systems), but it is an important engineering heuristic.

A one-step method applied to \( \dot{y}=\lambda y \) usually has the form \( y_{n+1}=R(z)y_n \) with \( z=h\lambda \), where \( R(z) \) is the stability function. The method is stable for a given \( z \) if \( |R(z)| \le 1 \).

Explicit Euler. From

\[ y_{n+1} = y_n + h\lambda y_n = (1+z)y_n, \]

we get the stability function

\[ R_{\mathrm{EE} }(z)=1+z. \]

For negative real \( \lambda \) (so \( z \in \mathbb{R}_{-} \)), stability requires

\[ |1+z| \le 1 \;\;\Longrightarrow\;\; -2 \le z \le 0. \]

Therefore for \( \lambda < 0 \), the step size must satisfy

\[ 0 < h \le \frac{2}{|\lambda|}. \]

If a model contains a fast stable mode such as \( \lambda=-1000 \), explicit Euler needs \( h \le 0.002 \), even if the slow physics evolves over seconds. This is the classical stiffness bottleneck.

Higher-order explicit Runge–Kutta methods improve accuracy, but their stability regions are still bounded. Hence they also suffer when \( h\lambda \) falls far into the left-half plane.

Implicit methods place the unknown state \( y_{n+1} \) inside the right-hand side, which leads to an algebraic equation at each step. The extra work per step is the price paid for much stronger stability properties.

4.1 Backward Euler (Implicit Euler)

Backward Euler is defined by

\[ y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}). \]

For the test equation, this becomes

\[ y_{n+1} = y_n + h\lambda y_{n+1} \;\;\Longrightarrow\;\; (1-z)y_{n+1}=y_n \;\;\Longrightarrow\;\; R_{\mathrm{BE} }(z)=\frac{1}{1-z}. \]

Proof sketch (A-stability of Backward Euler).

Let \( z=x+iy \) with \( x=\Re(z)\le 0 \). Then

\[ |1-z|^2 = |(1-x)-iy|^2 = (1-x)^2 + y^2 \ge 1. \]

Hence

\[ |R_{\mathrm{BE} }(z)| = \frac{1}{|1-z|} \le 1, \qquad \forall z \text{ with } \Re(z)\le 0. \]

Therefore Backward Euler is A-stable: its stability region contains the entire left-half complex plane.

It is also L-stable because

\[ \lim_{z \to -\infty} R_{\mathrm{BE} }(z)=0, \]

which means very fast decaying modes are strongly damped numerically.

4.2 Trapezoidal Rule (Implicit Midpoint Family Variant)

The trapezoidal rule (also called Crank–Nicolson in many contexts) is

\[ y_{n+1}=y_n+\frac{h}{2}\Big(f(t_n,y_n)+f(t_{n+1},y_{n+1})\Big). \]

For the test equation:

\[ y_{n+1}=y_n+\frac{z}{2}(y_n+y_{n+1}) \;\;\Longrightarrow\;\; \left(1-\frac{z}{2}\right)y_{n+1}=\left(1+\frac{z}{2}\right)y_n, \]

\[ R_{\mathrm{TR} }(z)=\frac{1+\frac{z}{2} }{1-\frac{z}{2} }. \]

Proof sketch (A-stability of Trapezoidal Rule).

Let \( z=x+iy \) with \( x\le 0 \). Then

\[ \left|1+\frac{z}{2}\right|^2 = \left(1+\frac{x}{2}\right)^2+\left(\frac{y}{2}\right)^2, \]

\[ \left|1-\frac{z}{2}\right|^2 = \left(1-\frac{x}{2}\right)^2+\left(\frac{y}{2}\right)^2. \]

Since \( x\le 0 \), we have \( \left(1+\frac{x}{2}\right)^2 \le \left(1-\frac{x}{2}\right)^2 \), so

\[ |R_{\mathrm{TR} }(z)| \le 1. \]

Thus trapezoidal is also A-stable. However, it is not L-stable because

\[ \lim_{z \to -\infty} R_{\mathrm{TR} }(z)=-1, \]

so very fast modes are not damped as aggressively as with Backward Euler.

4.3 BDF Family (High-Level)

Backward Differentiation Formulas (BDF) are multi-step implicit methods especially useful for stiff ODEs. For example, BDF2 is

\[ \frac{3y_{n+1}-4y_n+y_{n-1} }{2h}=f(t_{n+1},y_{n+1}). \]

Production solvers (e.g., BDF-based solvers) often use variable order and variable step size, combined with local error control and Jacobian reuse strategies.

The key computational task in implicit simulation is solving a nonlinear algebraic equation for \( y_{n+1} \). For Backward Euler:

\[ y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}). \]

Define the residual

\[ F(\mathbf{y}) = \mathbf{y} - \mathbf{y}_n - h\,\mathbf{f}(t_{n+1},\mathbf{y}). \]

Then we must solve \( F(\mathbf{y}_{n+1})=\mathbf{0} \). Newton’s method updates an iterate \( \mathbf{y}^{(m)} \) using

\[ \mathbf{J}_F(\mathbf{y}^{(m)})\,\Delta \mathbf{y}^{(m)} = -F(\mathbf{y}^{(m)}), \]

\[ \mathbf{y}^{(m+1)} = \mathbf{y}^{(m)} + \Delta \mathbf{y}^{(m)}, \]

where

\[ \mathbf{J}_F(\mathbf{y}) = \mathbf{I} - h\,\mathbf{J}_f(t_{n+1},\mathbf{y}), \qquad \mathbf{J}_f = \frac{\partial \mathbf{f} }{\partial \mathbf{y} }. \]

This is why Jacobians and linear solvers are central in stiff ODE integration. For large systems, sparse factorizations and Jacobian reuse can dominate performance.

flowchart TD
  A["Given y_n and step h"] --> B["Form residual F(y) = y - y_n - h f(t_{n+1}, y)"]
  B --> C["Initial guess for y_{n+1}"]
  C --> D["Newton iteration"]
  D --> E["Build J_F = I - h J_f"]
  E --> F["Solve linear system J_F * delta = -F"]
  F --> G["Update y = y + delta"]
  G --> H["Converged?"]
  H -->|no| D
  H -->|yes| I["Accept step and estimate error"]
        

Important engineering point: implicit methods are not “magic.” They convert step-size stability limits into algebraic solve effort. The practical win is that this trade-off is often favorable for stiff systems.

In Chapter 15 Lesson 2, you learned Runge–Kutta methods and step-size concepts. This lesson adds the next decision layer: when to switch to stiffness-aware methods.

Typical signs of stiffness in practice:

• Explicit solver takes extremely many tiny steps despite smooth state trajectories.
• Step size collapses after adding a highly damped subsystem or parasitic mode.
• A model with stable physics becomes numerically unstable unless \( h \) is tiny.
• Reaction/thermal/electrical transients happen on microseconds while control response is milliseconds/seconds.

Method-selection intuition (high level):

• Nonstiff: explicit RK methods (e.g., RK45) are efficient and easy.
• Mild-to-moderate stiffness: trapezoidal / implicit RK / Rosenbrock-type methods can work well.
• Strong stiffness: Backward Euler (robust but low-order) or variable-order BDF methods are common.

Libraries and tools commonly used in control/system dynamics simulation:

• Python: SciPy (solve_ivp with Radau, BDF)
• C++: SUNDIALS (CVODE), Boost.Odeint (nonstiff and some semi-implicit workflows), Eigen for linear algebra
• Java: Hipparchus / Apache Commons Math ODE packages, or custom implicit solvers for teaching
• MATLAB: ode15s, ode23s, ode23t, ode23tb
• Simulink: variable-step stiff solvers (especially ode15s) with tolerance control
• Wolfram Mathematica: NDSolve with stiffness-switching and BDF-based methods

This example compares RK45 with stiff-aware methods (Radau, BDF) on the Robertson stiff kinetics problem, then demonstrates Backward Euler on the scalar test equation.

Code file: Chapter15_Lesson3.py

"""
Chapter 15, Lesson 3 (System Dynamics)
Stiff systems and implicit methods: compare explicit and implicit solvers.
"""
import numpy as np
from scipy.integrate import solve_ivp
from time import perf_counter

def robertson(t, y):
    y1, y2, y3 = y
    return [
        -0.04*y1 + 1.0e4*y2*y3,
        0.04*y1 - 1.0e4*y2*y3 - 3.0e7*y2*y2,
        3.0e7*y2*y2
    ]

# Compare a nonstiff-oriented method vs stiff methods on a classic stiff IVP
t_span = (0.0, 1.0e2)
y0 = [1.0, 0.0, 0.0]

methods = ["RK45", "Radau", "BDF"]
for method in methods:
    t0 = perf_counter()
    sol = solve_ivp(robertson, t_span, y0, method=method, rtol=1e-6, atol=1e-10)
    dt = perf_counter() - t0
    print(f"Method={method:>5s} | success={sol.success} | nfev={sol.nfev:6d} | njev={getattr(sol,'njev',0):4d} | time={dt:.4f}s")
    print("  y(tf) =", sol.y[:, -1])

print("\n--- Backward Euler on scalar stiff test equation y' = lambda*y ---")
lam = -1000.0
h = 0.05
N = 10
y = 1.0
t = 0.0
for n in range(N):
    # For y' = lam*y, backward Euler gives y_{n+1} = y_n / (1 - h*lam)
    y = y / (1.0 - h*lam)
    t += h
    print(f"n={n+1:2d}, t={t:.3f}, y={y:.6e}")
      

In a stiff case, RK45 usually requires many more function evaluations than Radau or BDF. This directly demonstrates the “stability-limited step size” phenomenon discussed in Sections 2–4.

This teaching example implements Implicit Euler for a stiff 2-state linear system using a hand-written 2x2 linear solver. In production, one typically uses SUNDIALS (CVODE) and a robust linear algebra backend.

Code file: Chapter15_Lesson3.cpp

// Chapter15_Lesson3.cpp
// Implicit Euler for a stiff linear 2x2 system x' = A x
// Build example:
//   g++ -O2 -std=c++17 Chapter15_Lesson3.cpp -o stiff_demo
#include <iostream>
#include <array>
#include <iomanip>

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double,2>,2>;

Vec2 matVec(const Mat2& M, const Vec2& x) {
    return {
        M[0][0]*x[0] + M[0][1]*x[1],
        M[1][0]*x[0] + M[1][1]*x[1]
    };
}

// Solve 2x2 linear system M z = b
Vec2 solve2x2(const Mat2& M, const Vec2& b) {
    double a = M[0][0], c = M[0][1];
    double d = M[1][0], e = M[1][1];
    double det = a*e - c*d;
    if (std::abs(det) < 1e-15) {
        throw std::runtime_error("Singular matrix in solve2x2");
    }
    return {
        ( e*b[0] - c*b[1]) / det,
        (-d*b[0] + a*b[1]) / det
    };
}

int main() {
    // Stiff eigenvalues approximately -1 and -1000 (similarity-transformed diagonal)
    Mat2 A = { { {-1000.0, 999.0},
               {   0.0,   -1.0} } };

    double h = 0.02;
    int N = 50;
    Vec2 x{1.0, 1.0};

    // Implicit Euler: (I - h A) x_{n+1} = x_n
    Mat2 M = { { {1.0 - h*A[0][0], -h*A[0][1]},
               {-h*A[1][0],      1.0 - h*A[1][1]} } };

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Implicit Euler for stiff linear system\n";
    std::cout << "h = " << h << ", steps = " << N << "\n";
    double t = 0.0;
    for (int n = 0; n < N; ++n) {
        x = solve2x2(M, x);
        t += h;
        std::cout << "n=" << std::setw(3) << (n+1)
                  << "  t=" << std::setw(8) << t
                  << "  x=[" << std::setw(12) << x[0]
                  << ", " << std::setw(12) << x[1] << "]\n";
    }

    // For production-grade stiff solvers in C++, see SUNDIALS (CVODE) or odeint wrappers.
    return 0;
}
      

For linear systems, implicit Euler is especially simple because the step equation is linear. For nonlinear systems, the same pattern remains, but the matrix must be rebuilt (or approximated) inside a nonlinear iteration.

This example demonstrates the algebraic-solve idea directly in Java using Newton’s method for a scalar stiff nonlinear ODE. It is a good educational bridge before using larger libraries such as Hipparchus.

Code file: Chapter15_Lesson3.java

// Chapter15_Lesson3.java
// Backward Euler with Newton iterations for a scalar nonlinear stiff ODE:
//   y' = -1000 (y - cos t) - sin t
// Exact solution y(t) = cos t for consistent initial condition y(0)=1.

public class Chapter15_Lesson3 {

    static double f(double t, double y) {
        return -1000.0 * (y - Math.cos(t)) - Math.sin(t);
    }

    static double dfdY(double t, double y) {
        return -1000.0;
    }

    static double backwardEulerNewtonStep(double tNext, double yPrev, double h) {
        double y = yPrev; // initial guess
        int maxIter = 20;
        double tol = 1e-12;

        for (int k = 0; k < maxIter; k++) {
            double F = y - yPrev - h * f(tNext, y);
            double dF = 1.0 - h * dfdY(tNext, y);
            double delta = -F / dF;
            y += delta;
            if (Math.abs(delta) < tol) {
                break;
            }
        }
        return y;
    }

    public static void main(String[] args) {
        double h = 0.05;
        int N = 40;
        double t = 0.0;
        double y = 1.0;

        System.out.println("Backward Euler + Newton (Java)");
        System.out.printf("%5s %12s %18s %18s%n", "n", "t", "y_num", "y_exact");

        for (int n = 1; n <= N; n++) {
            double tNext = n * h;
            y = backwardEulerNewtonStep(tNext, y, h);
            double yExact = Math.cos(tNext);
            System.out.printf("%5d %12.6f %18.10f %18.10f%n", n, tNext, y, yExact);
            t = tNext;
        }

        // For production Java libraries, consider Hipparchus ODE solvers for broader method support.
    }
}
      

The derivative dfdY plays the role of the Jacobian in the scalar case. In vector systems, the same logic becomes a matrix solve at each Newton iteration.

MATLAB provides several stiff solvers. The most widely used general-purpose stiff solver is ode15s (variable-order, variable-step, based on numerical differentiation formulas). In Simulink, the same lesson applies: choose a stiff variable-step solver when fast stable modes dominate numerical behavior.

Code file: Chapter15_Lesson3.m

% Chapter15_Lesson3.m
% Stiff systems and implicit methods (MATLAB / Simulink-oriented demo)
% 1) Compare ode45 and ode15s on the Robertson problem
% 2) Show recommended Simulink solver settings for stiff models

function Chapter15_Lesson3()
    clc; clear;

    y0 = [1; 0; 0];
    tspan = [0 1e2];

    fprintf('--- Robertson problem: ode45 vs ode15s ---\n');
    opts = odeset('RelTol',1e-6,'AbsTol',1e-10);

    tic;
    [t45, y45] = ode45(@robertson, tspan, y0, opts);
    time45 = toc;

    tic;
    [t15s, y15s] = ode15s(@robertson, tspan, y0, opts);
    time15s = toc;

    fprintf('ode45 : steps = %d, elapsed = %.4f s\n', numel(t45), time45);
    fprintf('ode15s: steps = %d, elapsed = %.4f s\n', numel(t15s), time15s);
    fprintf('Final state ode45 : [%g, %g, %g]\n', y45(end,1), y45(end,2), y45(end,3));
    fprintf('Final state ode15s: [%g, %g, %g]\n', y15s(end,1), y15s(end,2), y15s(end,3));

    figure;
    semilogx(t15s + 1e-12, y15s(:,1), '-', t15s + 1e-12, y15s(:,2), '--', t15s + 1e-12, y15s(:,3), '-.');
    grid on;
    xlabel('t'); ylabel('States');
    title('Robertson problem solved with ode15s');
    legend('y1','y2','y3','Location','best');

    % Simulink note (programmatic setting):
    % If you already have a Simulink model named "StiffPlant", configure solver by:
    % set_param('StiffPlant', 'SolverType', 'Variable-step');
    % set_param('StiffPlant', 'Solver', 'ode15s');
    % set_param('StiffPlant', 'RelTol', '1e-6');
    % set_param('StiffPlant', 'AbsTol', '1e-8');
    % set_param('StiffPlant', 'MaxStep', 'auto');
end

function dydt = robertson(~, y)
    dydt = zeros(3,1);
    dydt(1) = -0.04*y(1) + 1e4*y(2)*y(3);
    dydt(2) =  0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2;
    dydt(3) =  3e7*y(2)^2;
end
      

Simulink configuration note: in stiff models, also inspect zero-crossing settings, tolerances, and any algebraic loop diagnostics. Solver choice is only one part of a stable simulation setup.

Mathematica’s NDSolve provides automatic stiffness detection/switching and stiff methods such as BDF. This is useful for quickly validating mathematical models before lower-level implementations.

Code file: Chapter15_Lesson3.nb

(* Chapter15_Lesson3.nb
   Wolfram Mathematica code (plain text notebook-ready content) *)

ClearAll["Global`*"];

robertson = {
  y1'[t] == -0.04 y1[t] + 10^4 y2[t] y3[t],
  y2'[t] == 0.04 y1[t] - 10^4 y2[t] y3[t] - 3*10^7 y2[t]^2,
  y3'[t] == 3*10^7 y2[t]^2,
  y1[0] == 1, y2[0] == 0, y3[0] == 0
};

solStiff = NDSolve[
  robertson,
  {y1, y2, y3},
  {t, 0, 100},
  Method -> {"StiffnessSwitching"},
  AccuracyGoal -> 10,
  PrecisionGoal -> 10
];

Print["Final state at t=100 (StiffnessSwitching):"];
Print[{y1[100], y2[100], y3[100]} /. solStiff // N];

solBDF = NDSolve[
  robertson,
  {y1, y2, y3},
  {t, 0, 100},
  Method -> {"BDF"},
  MaxStepFraction -> 1/50
];

Plot[
  Evaluate[{y1[t], y2[t], y3[t]} /. solBDF],
  {t, 0, 100},
  PlotRange -> All,
  PlotLegends -> {"y1", "y2", "y3"},
  GridLines -> Automatic,
  AxesLabel -> {"t", "state"},
  PlotLabel -> "Robertson Problem (BDF)"
];

(* Scalar stiff test equation with Backward Euler *)
lambda = -1000.0;
h = 0.05;
nSteps = 10;
y = 1.0;
Table[
  y = y/(1 - h*lambda);
  {n, n*h, y},
  {n, 1, nSteps}
] // TableForm
      

Problem 1 (Explicit Euler stability limit): Consider the stiff test equation \( \dot{y} = -500 y \). Find the maximum stable step size for Explicit Euler and explain why this can be restrictive.

Solution: Explicit Euler requires

\[ |1+h\lambda| \le 1. \]

With \( \lambda=-500 \), this gives

\[ -2 \le -500h \le 0 \;\;\Longrightarrow\;\; 0 \le h \le \frac{2}{500}=0.004. \]

So the step size must be at most \( 0.004 \), even if the slow behavior of the physical system evolves on a much larger time scale. This is the numerical signature of stiffness.

Problem 2 (Backward Euler on the test equation): Apply Backward Euler to \( \dot{y}=\lambda y \) and derive the amplification factor. Show why it is stable for all \( \Re(\lambda)\le 0 \).

Solution: Backward Euler gives

\[ y_{n+1}=y_n+h\lambda y_{n+1}. \]

Rearranging:

\[ (1-h\lambda)y_{n+1}=y_n, \qquad y_{n+1}=\frac{1}{1-h\lambda}y_n. \]

Thus

\[ R(z)=\frac{1}{1-z}, \qquad z=h\lambda. \]

If \( \Re(z)\le 0 \), then \( |1-z|\ge 1 \), so \( |R(z)|\le 1 \). Therefore Backward Euler is A-stable.

Problem 3 (Linear-system implicit step): For the linear system \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), derive the Implicit Euler update in matrix form.

Solution: Implicit Euler writes

\[ \mathbf{x}_{n+1}=\mathbf{x}_n+h\mathbf{A}\mathbf{x}_{n+1}. \]

Collecting terms on the left:

\[ (\mathbf{I}-h\mathbf{A})\mathbf{x}_{n+1}=\mathbf{x}_n. \]

Hence the step is obtained by a linear solve:

\[ \mathbf{x}_{n+1}=(\mathbf{I}-h\mathbf{A})^{-1}\mathbf{x}_n. \]

This formula explains the central role of matrix factorization in stiff linear simulation.

Problem 4 (Newton residual for Backward Euler): For a nonlinear ODE \( \dot{\mathbf{x} }=\mathbf{f}(t,\mathbf{x}) \), define the nonlinear residual for Backward Euler and derive its Jacobian.

Solution: Define

\[ F(\mathbf{x})=\mathbf{x}-\mathbf{x}_n-h\mathbf{f}(t_{n+1},\mathbf{x}). \]

The Jacobian of the residual is

\[ \mathbf{J}_F(\mathbf{x}) = \frac{\partial F}{\partial \mathbf{x} } = \mathbf{I}-h\frac{\partial \mathbf{f} }{\partial \mathbf{x} }(t_{n+1},\mathbf{x}). \]

Newton’s method then solves \( \mathbf{J}_F(\mathbf{x}^{(m)})\Delta \mathbf{x}^{(m)}=-F(\mathbf{x}^{(m)}) \) and updates \( \mathbf{x}^{(m+1)}=\mathbf{x}^{(m)}+\Delta \mathbf{x}^{(m)} \).

Problem 5 (A-stable but not L-stable): The trapezoidal rule has \( R(z)=\frac{1+z/2}{1-z/2} \). Show that it is not L-stable.

Solution: L-stability requires

\[ \lim_{z \to -\infty} R(z)=0. \]

For trapezoidal,

\[ \lim_{z \to -\infty}\frac{1+z/2}{1-z/2} = \lim_{z \to -\infty}\frac{z/2}{-z/2} = -1 \neq 0. \]

Therefore trapezoidal is not L-stable. It is still A-stable, but it does not strongly damp the fastest decaying modes.

Stiffness is a numerical phenomenon caused by fast stable modes that severely restrict explicit step sizes. The test-equation framework shows this through the stability function \( R(z) \). Implicit methods improve stability by solving algebraic equations each step: Backward Euler is A-stable and L-stable (very robust for stiff decay), while trapezoidal is A-stable but not L-stable. In practice, stiff-aware solvers (BDF/Radau/ode15s/StiffnessSwitching) are preferred when explicit methods become stability-limited.