function [A, B, B_w, x_star, u_star, Q_star] = pke_linearize(plant, x_star, u_star, Q_star)
% PKE_LINEARIZE  Numerical Jacobians of the PKE+T/H RHS at an operating point.
%
%  Computes (A, B, B_w) such that, for small perturbations dx, du, dw
%  around (x_star, u_star, Q_star):
%
%     dx/dt ~ A*dx + B*du + B_w*dw,   where  w = Q_sg.
%
%  Uses central finite differences.  Step sizes are scaled to each
%  variable's magnitude; fall back to an absolute epsilon near zero.
%
%  Inputs:
%    plant   - parameter struct from pke_params()
%    x_star  - (optional) 10x1 linearization point; default = operating x0
%    u_star  - (optional) scalar control; default = 0
%    Q_star  - (optional) scalar SG heat removal [W]; default = plant.P0
%
%  Outputs:
%    A       - 10x10 state Jacobian df/dx
%    B       - 10x1  input Jacobian df/du
%    B_w     - 10x1  disturbance Jacobian df/d(Q_sg)
%    x_star, u_star, Q_star - echoed back for the caller's records

    if nargin < 2 || isempty(x_star), x_star = pke_initial_conditions(plant); end
    if nargin < 3 || isempty(u_star), u_star = 0;                              end
    if nargin < 4 || isempty(Q_star), Q_star = plant.P0;                       end

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

    % Wrap the RHS as a function of (x, u, w) with plant fixed.
    f = @(x, u, w) pke_th_rhs(0, x, plant, @(t) w, u);

    f0 = f(x_star, u_star, Q_star);   %#ok<NASGU>  % sanity slot

    % --- A: df/dx via central differences, column by column ---
    A = zeros(n, n);
    for k = 1:n
        h  = max(eps_rel * abs(x_star(k)), eps_abs);
        xp = x_star;  xp(k) = xp(k) + h;
        xm = x_star;  xm(k) = xm(k) - h;
        A(:, k) = (f(xp, u_star, Q_star) - f(xm, u_star, Q_star)) / (2*h);
    end

    % --- B: df/du ---
    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)) / (2*h);

    % --- B_w: df/dQ_sg ---
    h   = max(eps_rel * abs(Q_star), 1.0);   % Q_star ~ 1e9, abs floor at 1 W
    B_w = (f(x_star, u_star, Q_star + h) - f(x_star, u_star, Q_star - h)) / (2*h);

end