"""
    pke_linearize(plant; x_star=nothing, u_star=0.0, Q_star=nothing)

Numerical Jacobians via central finite differences.  Returns
`(A, B, B_w, x_star, u_star, Q_star)` such that for small (dx, du, dw):

    dx/dt ≈ A dx + B du + B_w dw,   where w = Q_sg.

Defaults: `x_star` = operating-point steady state, `u_star = 0`,
`Q_star = P0`.
"""
function pke_linearize(plant; x_star=nothing, u_star=0.0, Q_star=nothing)
    x_star === nothing && (x_star = pke_initial_conditions(plant))
    Q_star === nothing && (Q_star = plant.P0)

    n = length(x_star)
    eps_rel = 1e-6
    eps_abs = 1e-8

    f = (x, u, w) -> pke_th_rhs(x, 0.0, plant, t -> w, u)

    A = zeros(n, n)
    for k in 1:n
        h = max(eps_rel * abs(x_star[k]), eps_abs)
        xp = copy(x_star);  xp[k] += h
        xm = copy(x_star);  xm[k] -= h
        A[:, k] = (f(xp, u_star, Q_star) - f(xm, u_star, Q_star)) ./ (2h)
    end

    h = max(eps_rel * abs(u_star), eps_abs)
    B = (f(x_star, u_star + h, Q_star) - f(x_star, u_star - h, Q_star)) ./ (2h)

    h = max(eps_rel * abs(Q_star), 1.0)
    B_w = (f(x_star, u_star, Q_star + h) - f(x_star, u_star, Q_star - h)) ./ (2h)

    return A, B, B_w, x_star, u_star, Q_star
end