function [t, X, U] = pke_solver(plant, Q_sg, ctrl_fn, ref, tspan)
% PKE_SOLVER  Solve the coupled PKE + T/H system in closed loop with a mode controller.
%
%  Inputs:
%    plant    - parameter struct from pke_params()
%    Q_sg     - function handle Q_sg(t) returning SG heat removal [W]
%    ctrl_fn  - function handle u = ctrl_fn(t, x, plant, ref)
%    ref      - struct of per-mode setpoints (e.g. ref.T_avg); [] if unused
%    tspan    - [t_start, t_end] in seconds
%
%  Outputs:
%    t  - time vector (M x 1) from the ODE solver
%    X  - state matrix (M x 10): columns are [n, C1..C6, T_f, T_c, T_cold]
%    U  - control input trajectory (M x 1): u evaluated at each (t(k), X(k,:))

    CtoF = @(T) T * 9/5 + 32;

    %% ===== Print Steady-State =====
    fprintf('=== Steady-State Conditions ===\n');
    fprintf('  beta       = %.5f  (%.1f pcm)\n', plant.beta, plant.beta*1e5);
    fprintf('  Lambda     = %.1e s\n', plant.Lambda);
    fprintf('  P0         = %.0f MWth\n', plant.P0/1e6);
    fprintf('  T_cold     = %.1f F\n', CtoF(plant.T_cold0));
    fprintf('  T_avg      = %.1f F\n', CtoF(plant.T_c0));
    fprintf('  T_hot      = %.1f F\n', CtoF(plant.T_hot0));
    fprintf('  T_fuel     = %.1f F\n', CtoF(plant.T_f0));
    fprintf('  Core dT    = %.1f F\n', plant.dT_core * 9/5);
    fprintf('================================\n\n');

    %% ===== Solve closed-loop =====
    x0      = pke_initial_conditions(plant);
    options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);

    rhs = @(t, x) pke_th_rhs(t, x, plant, Q_sg, ctrl_fn(t, x, plant, ref));
    [t, X] = ode15s(rhs, tspan, x0, options);

    %% ===== Reconstruct control trajectory =====
    M = length(t);
    U = zeros(M, 1);
    for k = 1:M
        U(k) = ctrl_fn(t(k), X(k,:).', plant, ref);
    end

    %% ===== Print Final State =====
    n = X(end, 1);  T_f = X(end, 8);  T_c = X(end, 9);  T_cold = X(end, 10);
    T_hot = 2*T_c - T_cold;
    rho_tot = U(end) ...
            + plant.alpha_f*(T_f - plant.T_f0) ...
            + plant.alpha_c*(T_c - plant.T_c0);

    fprintf('=== Final State (t = %.0f s) ===\n', t(end));
    fprintf('  Power  = %.1f MWth  (%.3f x nominal)\n', n*plant.P0/1e6, n);
    fprintf('  T_fuel = %.1f F\n', CtoF(T_f));
    fprintf('  T_hot  = %.1f F\n', CtoF(T_hot));
    fprintf('  T_avg  = %.1f F\n', CtoF(T_c));
    fprintf('  T_cold = %.1f F\n', CtoF(T_cold));
    fprintf('  u      = %.4f $  (external reactivity from controller)\n', U(end)/plant.beta);
    fprintf('  rho    = %.4f $  (total, incl. T-feedback)\n', rho_tot/plant.beta);

end