In Lesson 1, we converted an n-th order scalar ODE into a first-order vector system by defining a state as a stack of derivatives. In this lesson we specialize that construction into a structured state-space representation called the companion (phase-variable) form. Its system matrix has a distinctive “shift” structure (ones on the superdiagonal), and the final row contains the ODE coefficients.

We restrict attention to a scalar input \( u(t) \) entering the ODE without derivatives (a standard starting point in modern control modeling). The multi-input/multi-output extensions are treated in Chapter 5, Lesson 3.

flowchart TD
  A["Given: n-th order ODE with coefficients"] --> B["Normalize (make leading coefficient = 1)"]
  B --> C["Define states: x1=y, x2=dy/dt, ..., xn=d^(n-1)y/dt^(n-1)"]
  C --> D["Write chain: x1dot=x2, x2dot=x3, ..., x(n-1)dot=xn"]
  D --> E["Use ODE to express xndot as linear combo of x's + input u"]
  E --> F["Assemble matrices A, B, C, D (companion structure)"]
        

Consider the scalar LTI ODE of order \( n \ge 1 \):

\[ y^{(n)}(t) + a_{n-1}y^{(n-1)}(t) + \cdots + a_1 \dot{y}(t) + a_0 y(t) = b\,u(t), \quad a_i, b \in \mathbb{R}. \]

If the ODE is given with a non-unit leading coefficient, e.g. \( \alpha_n y^{(n)} + \cdots = b u \) with \( \alpha_n \neq 0 \), divide by \( \alpha_n \) to obtain the normalized (monic) form above. This step is essential because the canonical companion structure assumes the coefficient of \( y^{(n)} \) equals 1.

The model is fully specified by the input \( u(t) \) and the initial conditions: \( y(0), \dot{y}(0), \ldots, y^{(n-1)}(0) \).

Define the state vector as the output and its first \( n-1 \) derivatives:

\[ \mathbf{x}(t) := \begin{bmatrix} x_1(t)\\ x_2(t)\\ \vdots\\ x_n(t) \end{bmatrix} := \begin{bmatrix} y(t)\\ \dot{y}(t)\\ \vdots\\ y^{(n-1)}(t) \end{bmatrix}. \]

Then, by differentiation of the definitions, we immediately obtain the “shift” relations:

\[ \dot{x}_1(t)=x_2(t),\quad \dot{x}_2(t)=x_3(t),\quad \ldots,\quad \dot{x}_{n-1}(t)=x_n(t). \]

The only remaining component is \( \dot{x}_n(t)=y^{(n)}(t) \), which we get from the ODE by solving for \( y^{(n)} \).

flowchart LR
  U["u(t)"] --> SUM["sum with state terms"]
  X1["x1 = y"] --> SUM
  X2["x2 = dy/dt"] --> SUM
  XN["xn = d^(n-1)y/dt^(n-1)"] --> SUM
  SUM --> XND["xndot = -a0*x1 - a1*x2 - ... - a(n-1)*xn + b*u"]
  X1 --> X1D["x1dot = x2"]
  X2 --> X2D["x2dot = x3"]
  XN --> XNM1["x(n-1)dot = xn"]
        

Solve the ODE for \( y^{(n)}(t) \):

\[ y^{(n)}(t) = -a_{n-1}y^{(n-1)}(t) - \cdots - a_1\dot{y}(t) - a_0 y(t) + b\,u(t). \]

Substitute \( y^{(k)}(t)=x_{k+1}(t) \) for \( k=0,\ldots,n-1 \):

\[ \dot{x}_n(t) = -a_0 x_1(t) - a_1 x_2(t) - \cdots - a_{n-1} x_n(t) + b\,u(t). \]

Collecting all equations, we obtain the state-space form \( \dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}u(t) \) with companion structure:

\[ \mathbf{A} := \begin{bmatrix} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_0 & -a_1 & -a_2 & \cdots & -a_{n-1} \end{bmatrix}, \quad \mathbf{B} := \begin{bmatrix} 0\\ 0\\ \vdots\\ 0\\ b \end{bmatrix}. \]

If we choose the output as the original variable \( y(t) \), then \( y(t)=x_1(t) \), hence:

\[ \mathbf{C} := \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}, \qquad \mathbf{D} := \begin{bmatrix} 0 \end{bmatrix}. \]

This realization is frequently called the phase-variable form because the state variables resemble “generalized coordinates” along a derivative chain.

We now prove that the companion-form state-space model is equivalent to the original ODE in the sense of producing the same output trajectory for the same input and consistent initial data.

Theorem (ODE ↔ companion state-space equivalence).

Let \( y(t) \) satisfy the monic ODE \( y^{(n)} + a_{n-1}y^{(n-1)} + \cdots + a_0 y = b u \). Define \( \mathbf{x}(t) = [y, \dot{y}, \ldots, y^{(n-1)}]^T \). Then \( \mathbf{x}(t) \) satisfies \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \) with companion \( \mathbf{A},\mathbf{B} \) from Section 4, and the output \( y=\mathbf{C}\mathbf{x} \) equals the original \( y(t) \). Conversely, if \( \mathbf{x}(t) \) satisfies the companion state equation and \( y(t)=x_1(t) \), then \( y(t) \) satisfies the ODE.

Proof.

(Forward direction.) By definition, \( x_1=y \) and \( x_{k+1}=y^{(k)} \) for \( k=1,\ldots,n-1 \). Differentiating gives \( \dot{x}_k = x_{k+1} \) for \( k=1,\ldots,n-1 \). Also \( \dot{x}_n = y^{(n)} \). Using the ODE solved for \( y^{(n)} \), we get \( \dot{x}_n = -a_0x_1 - a_1x_2 - \cdots - a_{n-1}x_n + b u \). These are exactly the components of \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \) with the stated companion matrices. Finally \( y=x_1=\mathbf{C}\mathbf{x} \).

(Reverse direction.) Assume \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \) and set \( y=x_1 \). The first \( n-1 \) equations imply \( \dot{x}_1=x_2, \dot{x}_2=x_3, \ldots, \dot{x}_{n-1}=x_n \), so by repeated differentiation \( x_{k+1}=y^{(k)} \) for \( k=1,\ldots,n-1 \). The last equation gives \( y^{(n)}=\dot{x}_n = -a_0x_1 - \cdots - a_{n-1}x_n + b u \). Substituting \( x_1=y, x_2=\dot{y}, \ldots, x_n=y^{(n-1)} \) yields the original ODE.

Therefore the two descriptions are equivalent. □

Initial-condition consistency. The state initial condition must satisfy: \( \mathbf{x}(0)= [y(0), \dot{y}(0), \ldots, y^{(n-1)}(0)]^T \).

Consider:

\[ \ddot{y}(t) + 3\dot{y}(t) + 2y(t) = u(t). \]

Here \( n=2 \), \( a_1=3 \), \( a_0=2 \), \( b=1 \). Define \( x_1=y \), \( x_2=\dot{y} \). Then:

\[ \dot{x}_1 = x_2,\qquad \dot{x}_2 = -2x_1 - 3x_2 + u. \]

Therefore:

\[ \mathbf{A}= \begin{bmatrix} 0 & 1\\ -2 & -3 \end{bmatrix}, \quad \mathbf{B}= \begin{bmatrix} 0\\ 1 \end{bmatrix}, \quad \mathbf{C}= \begin{bmatrix} 1 & 0 \end{bmatrix}, \quad \mathbf{D}= \begin{bmatrix} 0 \end{bmatrix}. \]

This is the canonical “chain-of-integrators + feedback” structure that will recur throughout modern control.

Recommended libraries for modern control workflows in Python: \( \) numpy, scipy (simulation/linear algebra), and optionally control (state-space utilities and canonical realizations). Below we build \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) directly from ODE coefficients and simulate the state equation with solve_ivp.

import numpy as np
from scipy.integrate import solve_ivp

def companion_from_ode(a, b=1.0):
    """
    Build companion (phase-variable) matrices for:
        y^(n) + a[n-1] y^(n-1) + ... + a[1] y' + a[0] y = b u
    Input:
        a: array-like length n, ordered [a0, a1, ..., a(n-1)]
        b: scalar
    Returns:
        A (n,n), B (n,1), C (1,n), D (1,1)
    """
    a = np.asarray(a, dtype=float)
    n = a.size
    A = np.zeros((n, n))
    # superdiagonal ones
    for i in range(n - 1):
        A[i, i + 1] = 1.0
    # last row: [-a0, -a1, ..., -a(n-1)]
    A[-1, :] = -a
    B = np.zeros((n, 1))
    B[-1, 0] = float(b)
    C = np.zeros((1, n))
    C[0, 0] = 1.0
    D = np.zeros((1, 1))
    return A, B, C, D

# Example: y'' + 3 y' + 2 y = u
A, B, C, D = companion_from_ode(a=[2.0, 3.0], b=1.0)

def u_of_t(t):
    # Example input: unit step
    return 1.0 if t >= 0.0 else 0.0

def f(t, x):
    x = np.asarray(x).reshape(-1, 1)
    dx = A @ x + B * u_of_t(t)
    return dx.flatten()

x0 = np.array([0.0, 0.0])  # y(0), y'(0)
t_span = (0.0, 5.0)
t_eval = np.linspace(t_span[0], t_span[1], 400)

sol = solve_ivp(f, t_span, x0, t_eval=t_eval, rtol=1e-9, atol=1e-12)
y = (C @ sol.y).flatten()

print("A=\n", A)
print("B=\n", B)
print("y(t) first/last:", y[0], y[-1])

# Optional: if python-control is installed
# from control import ss
# sys = ss(A, B, C, D)

The function companion_from_ode is suitable for arbitrary order \( n \). In later chapters, the same structure will be leveraged for design and analysis routines.

Common C++ choices for modern control numerics include: Eigen (linear algebra) and either a custom integrator (e.g., RK4) or boost::numeric::odeint for ODE integration. Below is a compact Eigen-based RK4 simulation.

#include <Eigen/Dense>
#include <iostream>
#include <vector>
#include <functional>

struct SS {
  Eigen::MatrixXd A;
  Eigen::VectorXd B;
  Eigen::RowVectorXd C;
  double D;
};

SS companion_from_ode(const std::vector<double>& a, double b) {
  const int n = static_cast<int>(a.size());
  SS sys;
  sys.A = Eigen::MatrixXd::Zero(n, n);
  for (int i = 0; i < n - 1; ++i) sys.A(i, i + 1) = 1.0;
  for (int j = 0; j < n; ++j) sys.A(n - 1, j) = -a[j];
  sys.B = Eigen::VectorXd::Zero(n);
  sys.B(n - 1) = b;
  sys.C = Eigen::RowVectorXd::Zero(n);
  sys.C(0) = 1.0;
  sys.D = 0.0;
  return sys;
}

Eigen::VectorXd rk4_step(const std::function<Eigen::VectorXd(double, const Eigen::VectorXd&)>& f,
                         double t, const Eigen::VectorXd& x, double h) {
  Eigen::VectorXd k1 = f(t, x);
  Eigen::VectorXd k2 = f(t + 0.5*h, x + 0.5*h*k1);
  Eigen::VectorXd k3 = f(t + 0.5*h, x + 0.5*h*k2);
  Eigen::VectorXd k4 = f(t + h, x + h*k3);
  return x + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
}

int main() {
  // Example: y'' + 3 y' + 2 y = u
  SS sys = companion_from_ode({2.0, 3.0}, 1.0);

  auto u_of_t = [](double t) -> double { return (t >= 0.0) ? 1.0 : 0.0; };

  auto f = [&](double t, const Eigen::VectorXd& x) -> Eigen::VectorXd {
    return sys.A * x + sys.B * u_of_t(t);
  };

  double t0 = 0.0, tf = 5.0, h = 1e-3;
  int N = static_cast<int>((tf - t0) / h);
  Eigen::VectorXd x(2);
  x << 0.0, 0.0; // [y(0), y'(0)]
  double t = t0;

  for (int k = 0; k < N; ++k) {
    x = rk4_step(f, t, x, h);
    t += h;
  }

  double y = sys.C * x + sys.D * u_of_t(t);
  std::cout << "Final y(tf) = " << y << std::endl;
  return 0;
}

This code directly mirrors the mathematical construction: build companion matrices, then integrate \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \). In later lessons, more sophisticated solvers and discretization choices will be discussed.

In Java, typical libraries are: EJML for matrix operations and Apache Commons Math for numerical integration. Below is a self-contained RK4 simulation using EJML matrices.

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;

import java.util.function.DoubleFunction;

public class CompanionSS {

  static class SS {
    DMatrixRMaj A;   // n x n
    DMatrixRMaj B;   // n x 1
    DMatrixRMaj C;   // 1 x n
    double D;        // scalar
  }

  static SS companionFromOde(double[] a, double b) {
    int n = a.length;
    SS sys = new SS();
    sys.A = new DMatrixRMaj(n, n);
    sys.B = new DMatrixRMaj(n, 1);
    sys.C = new DMatrixRMaj(1, n);
    sys.D = 0.0;

    // superdiagonal ones
    for (int i = 0; i < n - 1; i++) sys.A.set(i, i + 1, 1.0);

    // last row = -a
    for (int j = 0; j < n; j++) sys.A.set(n - 1, j, -a[j]);

    sys.B.set(n - 1, 0, b);
    sys.C.set(0, 0, 1.0);
    return sys;
  }

  static DMatrixRMaj f(SS sys, double t, DMatrixRMaj x, DoubleFunction<Double> uOfT) {
    DMatrixRMaj Ax = new DMatrixRMaj(x.numRows, 1);
    CommonOps_DDRM.mult(sys.A, x, Ax);

    DMatrixRMaj Bu = new DMatrixRMaj(x.numRows, 1);
    CommonOps_DDRM.scale(uOfT.apply(t), sys.B, Bu);

    CommonOps_DDRM.addEquals(Ax, Bu);
    return Ax; // dx/dt
  }

  static DMatrixRMaj rk4Step(SS sys, double t, DMatrixRMaj x, double h, DoubleFunction<Double> uOfT) {
    DMatrixRMaj k1 = f(sys, t, x, uOfT);

    DMatrixRMaj x2 = x.copy();
    CommonOps_DDRM.addEquals(x2, 0.5*h, k1);
    DMatrixRMaj k2 = f(sys, t + 0.5*h, x2, uOfT);

    DMatrixRMaj x3 = x.copy();
    CommonOps_DDRM.addEquals(x3, 0.5*h, k2);
    DMatrixRMaj k3 = f(sys, t + 0.5*h, x3, uOfT);

    DMatrixRMaj x4 = x.copy();
    CommonOps_DDRM.addEquals(x4, h, k3);
    DMatrixRMaj k4 = f(sys, t + h, x4, uOfT);

    DMatrixRMaj out = x.copy();
    // out = x + (h/6)*(k1 + 2k2 + 2k3 + k4)
    CommonOps_DDRM.addEquals(out, h/6.0, k1);
    CommonOps_DDRM.addEquals(out, h/3.0, k2);
    CommonOps_DDRM.addEquals(out, h/3.0, k3);
    CommonOps_DDRM.addEquals(out, h/6.0, k4);
    return out;
  }

  public static void main(String[] args) {
    // Example: y'' + 3 y' + 2 y = u
    SS sys = companionFromOde(new double[]{2.0, 3.0}, 1.0);

    DoubleFunction<Double> uOfT = (double t) -> (t >= 0.0) ? 1.0 : 0.0;

    DMatrixRMaj x = new DMatrixRMaj(2, 1);
    x.set(0, 0, 0.0); // y(0)
    x.set(1, 0, 0.0); // y'(0)

    double t = 0.0, tf = 5.0, h = 1e-3;
    int N = (int)Math.round((tf - t)/h);

    for (int k = 0; k < N; k++) {
      x = rk4Step(sys, t, x, h, uOfT);
      t += h;
    }

    // y = C x + D u
    DMatrixRMaj y = new DMatrixRMaj(1, 1);
    CommonOps_DDRM.mult(sys.C, x, y);
    y.set(0, 0, y.get(0, 0) + sys.D * uOfT.apply(t));
    System.out.println("Final y(tf) = " + y.get(0, 0));
  }
}

The construction is identical across languages; only the matrix and integration tooling changes.

MATLAB (with Control System Toolbox) supports state-space objects, simulation (lsim), and programmatic model construction. Here we construct the companion matrices directly from the ODE coefficients, then simulate the response to a step input.

function [A,B,C,D] = companion_from_ode(a,b)
% Build companion matrices for:
%   y^(n) + a(n) y^(n-1) + ... + a(2) y' + a(1) y = b u
% Input a is [a0 a1 ... a(n-1)] in the lesson's notation.

a = a(:).';              % row
n = length(a);
A = zeros(n,n);
A(1:n-1,2:n) = eye(n-1);
A(n,:) = -a;             % [-a0 ... -a(n-1)]
B = zeros(n,1); B(n) = b;
C = zeros(1,n); C(1) = 1;
D = 0;
end

% Example: y'' + 3 y' + 2 y = u
[A,B,C,D] = companion_from_ode([2 3], 1);

sys = ss(A,B,C,D);

t = linspace(0,5,400);
u = ones(size(t));       % step
x0 = [0;0];              % [y(0); y'(0)]
y = lsim(sys,u,t,x0);

disp(A); disp(B);
fprintf('y(tf) = %.6f\n', y(end));

Simulink realization idea (integrator chain).

Build a chain of n integrators whose outputs are \( x_n, x_{n-1}, \ldots, x_1 \) (or in forward order), then form the last derivative equation: \( \dot{x}_n = -a_0 x_1 - \cdots - a_{n-1} x_n + b u \) using a Sum block and Gain blocks. This exactly mirrors the companion structure derived above.

% Sketch: programmatic Simulink construction (outline)
% new_system('companion_model'); open_system('companion_model');
% Add blocks: In1 (u), Sum, Gains for -a_i, Integrators (n), Out1 (y)
% Wire:
%   x1dot = x2, x2dot = x3, ..., x(n-1)dot = xn
%   xndot = sum(-a_i * xi) + b*u
% Output y = x1
% Then set solver options and simulate via sim('companion_model')

Mathematica can symbolically form companion matrices and simulate the state equation using NDSolve. The following code constructs \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) and simulates a step input for the second-order example.

companionFromOde[a_List, b_: 1] := Module[{n, A, B, C, D},
  n = Length[a];
  A = ConstantArray[0, {n, n}];
  Do[A[[i, i + 1]] = 1, {i, 1, n - 1}];
  A[[n, All]] = -a;
  B = ConstantArray[0, {n, 1}];
  B[[n, 1]] = b;
  C = ConstantArray[0, {1, n}];
  C[[1, 1]] = 1;
  D = {{0}};
  {A, B, C, D}
];

(* Example: y'' + 3 y' + 2 y = u *)
{A, B, C, D} = companionFromOde[{2, 3}, 1];

u[t_] := 1;  (* step input *)
x10 = 0; x20 = 0;

sol = NDSolve[
  {
    x1'[t] == x2[t],
    x2'[t] == -2 x1[t] - 3 x2[t] + u[t],
    x1[0] == x10,
    x2[0] == x20
  },
  {x1, x2},
  {t, 0, 5}
];

y[t_] := x1[t] /. sol[[1]];
{A, B, C, D}
y[5]

For higher order \( n \), extend the state definitions and equations automatically from the constructed matrices, or use StateSpaceModel with \( \{\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\} \).

The following problems reinforce the derivation and ensure fluency in converting higher-order ODEs into the companion state-space form.

Problem 1 (Third-order construction). Consider the ODE: \( y^{(3)}(t) + 4y^{(2)}(t) + 5\dot{y}(t) + 2y(t) = 3u(t) \). Construct the companion matrices \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) with \( y=x_1 \).

Solution. Here \( n=3 \), and coefficients are: \( a_2=4 \), \( a_1=5 \), \( a_0=2 \), \( b=3 \). Define \( x_1=y \), \( x_2=\dot{y} \), \( x_3=y^{(2)} \).

\[ \dot{x}_1=x_2,\quad \dot{x}_2=x_3,\quad \dot{x}_3=-2x_1-5x_2-4x_3+3u. \]

\[ \mathbf{A}= \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -2 & -5 & -4 \end{bmatrix}, \quad \mathbf{B}= \begin{bmatrix} 0\\ 0\\ 3 \end{bmatrix}, \quad \mathbf{C}= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}, \quad \mathbf{D}= \begin{bmatrix} 0 \end{bmatrix}. \]

Problem 2 (Normalization to monic form). Given: \( 2y^{(4)}(t) + y^{(3)}(t) - 6y^{(2)}(t) + y(t) = u(t) \), (i) rewrite it in monic form, then (ii) write the last row of the companion matrix \( \mathbf{A} \) (you do not need to write the full \( \mathbf{A} \)).

Solution. Divide by 2:

\[ y^{(4)}(t) + \tfrac{1}{2}y^{(3)}(t) - 3y^{(2)}(t) + 0\cdot \dot{y}(t) + \tfrac{1}{2}y(t) = \tfrac{1}{2}u(t). \]

Therefore \( a_3=\tfrac{1}{2} \), \( a_2=-3 \), \( a_1=0 \), \( a_0=\tfrac{1}{2} \). The last row is:

\[ \mathbf{A}_{4,\cdot} = \begin{bmatrix} -a_0 & -a_1 & -a_2 & -a_3 \end{bmatrix} = \begin{bmatrix} -\tfrac{1}{2} & 0 & 3 & -\tfrac{1}{2} \end{bmatrix}. \]

Problem 3 (Recovering the ODE from a companion form). Suppose a state equation is given with: \( \mathbf{A}= \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -7 & 2 & -5 \end{bmatrix} \), \( \mathbf{B}= \begin{bmatrix} 0\\ 0\\ 4 \end{bmatrix} \), and \( y=x_1 \). Write the corresponding scalar ODE.

Solution. In companion form, the last row equals \( [-a_0, -a_1, -a_2] \). Hence: \( -a_0=-7 \Rightarrow a_0=7 \), \( -a_1=2 \Rightarrow a_1=-2 \), \( -a_2=-5 \Rightarrow a_2=5 \), and \( b=4 \).

\[ y^{(3)}(t) + 5y^{(2)}(t) - 2\dot{y}(t) + 7y(t) = 4u(t). \]

Problem 4 (Initial condition mapping). For the ODE \( y^{(3)} + a_2 y^{(2)} + a_1 \dot{y} + a_0 y = b u \), prove that the unique state initial condition consistent with ODE initial data \( y(0),\dot{y}(0),y^{(2)}(0) \) is \( \mathbf{x}(0)=[y(0),\dot{y}(0),y^{(2)}(0)]^T \).

Solution. By definition in phase variables, \( x_1(t)=y(t) \), \( x_2(t)=\dot{y}(t) \), \( x_3(t)=y^{(2)}(t) \). Evaluating at \( t=0 \) gives \( x_1(0)=y(0) \), \( x_2(0)=\dot{y}(0) \), \( x_3(0)=y^{(2)}(0) \). Any other initial state would contradict at least one of these defining identities, hence is inconsistent.

Problem 5 (General n, compact derivation). Let \( \mathbf{x}=[y,\dot{y},\ldots,y^{(n-1)}]^T \). Show that \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \) with the companion \( \mathbf{A},\mathbf{B} \) derived in Section 4, by writing the component-wise relations and collecting them into matrix form.

Solution. For \( i=1,\ldots,n-1 \), definitions imply \( \dot{x}_i=x_{i+1} \), which forms the superdiagonal ones in \( \mathbf{A} \). The ODE implies \( \dot{x}_n=y^{(n)}=-a_0x_1-\cdots-a_{n-1}x_n+b u \), which forms the last row of \( \mathbf{A} \) and the last entry of \( \mathbf{B} \). Writing the \( n \) scalar equations as a single vector equation yields the stated matrices.

We derived the companion (phase-variable) state-space form directly from a monic n-th order scalar LTI ODE. The key idea is to define the state as the output and its derivatives, producing a “shift” structure in \( \mathbf{A} \) and placing the ODE coefficients in the last row. We proved bidirectional equivalence between the ODE and the state-space model and implemented the construction across major computational environments used in modern control.