In multivariable dynamics, we often encounter coupled differential equations where the derivative of one variable depends on many others. Writing these equations component-by-component is correct but cumbersome. The vector-matrix form:

\[ \dot{\mathbf{x}}(t) = \mathbf{A}(t)\mathbf{x}(t) + \mathbf{b}(t), \quad \mathbf{x}(t)\in\mathbb{R}^n,\; \mathbf{A}(t)\in\mathbb{R}^{n\times n},\; \mathbf{b}(t)\in\mathbb{R}^n \]

makes structure explicit: coupling is encoded in \( \mathbf{A}(t) \), and forcing is encoded in \( \mathbf{b}(t) \). Later, this same form will become the standard representation for internal system dynamics.

flowchart TD
  S["Given coupled scalar ODEs"] --> M["Stack variables into vector x(t)"]
  M --> A["Collect coefficients into matrix A(t)"]
  A --> B["Collect forcing into vector b(t)"]
  B --> V["Obtain compact form: x_dot = A(t) x + b(t)"]
        

We distinguish two common cases:

• \( \mathbf{A}(t) \) varies with time (linear time-varying).
• \( \mathbf{A} \) is constant (linear time-invariant), which will be central in subsequent chapters.

Let \( \mathbf{x}:[t_0,t_1]\rightarrow\mathbb{R}^n \) be differentiable with components \( x_i(t) \). The derivative is defined componentwise:

\[ \mathbf{x}(t)=\begin{bmatrix}x_1(t)\\ \vdots\\ x_n(t)\end{bmatrix} \quad\Longrightarrow\quad \dot{\mathbf{x}}(t)=\frac{d\mathbf{x}}{dt}(t)=\begin{bmatrix}\dot{x}_1(t)\\ \vdots\\ \dot{x}_n(t)\end{bmatrix}. \]

Similarly, for a matrix-valued function \( \mathbf{A}(t)=[a_{ij}(t)] \), define:

\[ \dot{\mathbf{A}}(t)=\frac{d\mathbf{A}}{dt}(t)=\left[\dot{a}_{ij}(t)\right]. \]

Linearity of differentiation. For scalars \( \alpha,\beta\in\mathbb{R} \) and differentiable vectors \( \mathbf{x}(t),\mathbf{y}(t) \):

\[ \frac{d}{dt}\left(\alpha\mathbf{x}(t)+\beta\mathbf{y}(t)\right) =\alpha\dot{\mathbf{x}}(t)+\beta\dot{\mathbf{y}}(t). \]

Product rule (matrix-vector). If \( \mathbf{A}(t) \) and \( \mathbf{x}(t) \) are differentiable, then:

\[ \frac{d}{dt}\left(\mathbf{A}(t)\mathbf{x}(t)\right)=\dot{\mathbf{A}}(t)\mathbf{x}(t)+\mathbf{A}(t)\dot{\mathbf{x}}(t). \]

Proof (componentwise). The \( i \)-th component of \( \mathbf{A}(t)\mathbf{x}(t) \) is \( \sum_{j=1}^n a_{ij}(t)x_j(t) \). Differentiate using the scalar product rule:

\[ \frac{d}{dt}\sum_{j=1}^n a_{ij}(t)x_j(t) =\sum_{j=1}^n \dot{a}_{ij}(t)x_j(t)+\sum_{j=1}^n a_{ij}(t)\dot{x}_j(t), \]

which is exactly the \( i \)-th component of \( \dot{\mathbf{A}}(t)\mathbf{x}(t)+\mathbf{A}(t)\dot{\mathbf{x}}(t) \).

Consider \( n \) first-order linear ODEs:

\[ \dot{x}_i(t) = \sum_{j=1}^n a_{ij}(t)x_j(t) + b_i(t),\quad i=1,\dots,n. \]

Define \( \mathbf{x}(t)=[x_1(t),\dots,x_n(t)]^\top \), \( \mathbf{A}(t)=[a_{ij}(t)] \), and \( \mathbf{b}(t)=[b_1(t),\dots,b_n(t)]^\top \). Then the system becomes:

\[ \dot{\mathbf{x}}(t) = \mathbf{A}(t)\mathbf{x}(t) + \mathbf{b}(t). \]

Equivalence claim. The vector equation holds if and only if all scalar equations hold.

Proof. Expanding \( \mathbf{A}(t)\mathbf{x}(t) \) gives \( [\sum_{j}a_{1j}(t)x_j(t),\dots,\sum_{j}a_{nj}(t)x_j(t)]^\top \). Adding \( \mathbf{b}(t) \) yields the right-hand sides of each scalar equation.

This reformulation is not cosmetic: it enables the use of matrix algebra (rank, null space, similarity transforms) to analyze structure and, later, to derive closed-form solutions systematically.

Define the (time-varying) linear differential operator \( \mathcal{L} \) acting on differentiable vector functions:

\[ \mathcal{L}[\mathbf{x}](t) \triangleq \dot{\mathbf{x}}(t) - \mathbf{A}(t)\mathbf{x}(t). \]

Then the forced system is \( \mathcal{L}[\mathbf{x}](t)=\mathbf{b}(t) \).

Proposition (linearity of \( \mathcal{L} \)). For scalars \( \alpha,\beta \) and differentiable \( \mathbf{x},\mathbf{y} \):

\[ \mathcal{L}[\alpha\mathbf{x}+\beta\mathbf{y}](t)=\alpha\mathcal{L}[\mathbf{x}](t)+\beta\mathcal{L}[\mathbf{y}](t). \]

Proof. By linearity of differentiation and matrix multiplication:

\[ \mathcal{L}[\alpha\mathbf{x}+\beta\mathbf{y}] = \alpha\dot{\mathbf{x}}+\beta\dot{\mathbf{y}}-\mathbf{A}(\alpha\mathbf{x}+\beta\mathbf{y}) = \alpha(\dot{\mathbf{x}}-\mathbf{A}\mathbf{x})+\beta(\dot{\mathbf{y}}-\mathbf{A}\mathbf{y}). \]

Corollary (superposition for homogeneous systems). If \( \mathbf{b}(t)=\mathbf{0} \) and \( \mathbf{x}_1,\mathbf{x}_2 \) solve \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \), then \( \alpha\mathbf{x}_1+\beta\mathbf{x}_2 \) is also a solution.

Corollary (affine structure for forced systems). If \( \mathbf{x}_p \) is one particular solution of \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \), then every solution has the form:

\[ \mathbf{x}(t)=\mathbf{x}_p(t)+\mathbf{x}_h(t), \quad \text{where } \dot{\mathbf{x}}_h(t)=\mathbf{A}(t)\mathbf{x}_h(t). \]

Proof. Let \( \mathbf{x} \) be any solution and define \( \mathbf{x}_h=\mathbf{x}-\mathbf{x}_p \). Then \( \dot{\mathbf{x}}_h=\dot{\mathbf{x}}-\dot{\mathbf{x}}_p=\mathbf{A}\mathbf{x}+\mathbf{b}-\left(\mathbf{A}\mathbf{x}_p+\mathbf{b}\right)=\mathbf{A}\mathbf{x}_h \).

Even without solving explicitly, we can rewrite the system as an integral equation. Integrate both sides from \( t_0 \) to \( t \):

\[ \int_{t_0}^{t}\dot{\mathbf{x}}(s)\,ds = \int_{t_0}^{t}\left(\mathbf{A}(s)\mathbf{x}(s)+\mathbf{b}(s)\right)\,ds. \]

By the fundamental theorem of calculus, \( \int_{t_0}^{t}\dot{\mathbf{x}}(s)\,ds=\mathbf{x}(t)-\mathbf{x}(t_0) \), so:

\[ \mathbf{x}(t)=\mathbf{x}(t_0)+\int_{t_0}^{t}\mathbf{A}(s)\mathbf{x}(s)\,ds + \int_{t_0}^{t}\mathbf{b}(s)\,ds. \]

This representation is the starting point for existence/uniqueness analysis (formalized later) and for constructing solution operators (developed next via matrix exponential and fundamental matrices).

From Chapter 2, you know that a change of coordinates can simplify matrix structure. Let \( \mathbf{T}\in\mathbb{R}^{n\times n} \) be invertible and define a new variable: \( \mathbf{z}(t)=\mathbf{T}^{-1}\mathbf{x}(t) \), i.e., \( \mathbf{x}(t)=\mathbf{T}\mathbf{z}(t) \).

Case A: constant transformation \( \mathbf{T} \). Differentiate \( \mathbf{x}(t)=\mathbf{T}\mathbf{z}(t) \):

\[ \dot{\mathbf{x}}(t)=\mathbf{T}\dot{\mathbf{z}}(t). \]

Substitute into \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \):

\[ \mathbf{T}\dot{\mathbf{z}}(t)=\mathbf{A}(t)\mathbf{T}\mathbf{z}(t)+\mathbf{b}(t) \quad\Longrightarrow\quad \dot{\mathbf{z}}(t)=\underbrace{\mathbf{T}^{-1}\mathbf{A}(t)\mathbf{T}}_{\mathbf{A}_z(t)}\mathbf{z}(t)+\underbrace{\mathbf{T}^{-1}\mathbf{b}(t)}_{\mathbf{b}_z(t)}. \]

Therefore, under a constant coordinate change, the system matrix transforms by similarity: \( \mathbf{A}_z(t)=\mathbf{T}^{-1}\mathbf{A}(t)\mathbf{T} \).

Case B: time-varying transformation \( \mathbf{T}(t) \) (advanced note). If \( \mathbf{x}(t)=\mathbf{T}(t)\mathbf{z}(t) \), then by the product rule:

\[ \dot{\mathbf{x}}(t)=\dot{\mathbf{T}}(t)\mathbf{z}(t)+\mathbf{T}(t)\dot{\mathbf{z}}(t), \]

and the transformed dynamics become:

\[ \dot{\mathbf{z}}(t)=\mathbf{T}(t)^{-1}\mathbf{A}(t)\mathbf{T}(t)\mathbf{z}(t) - \mathbf{T}(t)^{-1}\dot{\mathbf{T}}(t)\mathbf{z}(t) + \mathbf{T}(t)^{-1}\mathbf{b}(t). \]

The extra term \( -\mathbf{T}^{-1}\dot{\mathbf{T}} \) is purely from time-varying coordinates. We will mainly use constant \( \mathbf{T} \) transformations in the near term.

Closed-form solutions are not always convenient, and simulation is essential. Let \( t_k=t_0+kh \) with step size \( h>0 \) and approximate \( \mathbf{x}(t_k)\approx \mathbf{x}_k \).

Forward Euler (first order). Using \( \dot{\mathbf{x}}(t_k)\approx (\mathbf{x}_{k+1}-\mathbf{x}_k)/h \):

\[ \mathbf{x}_{k+1}=\mathbf{x}_k + h\left(\mathbf{A}(t_k)\mathbf{x}_k+\mathbf{b}(t_k)\right). \]

The local truncation error is \( O(h^2) \) (hence global error \( O(h) \)) under standard smoothness assumptions.

RK4 (fourth order). Define the vector field: \( \mathbf{f}(t,\mathbf{x})=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \). Then:

\[ \begin{aligned} \mathbf{k}_1 & = \mathbf{f}(t_k,\mathbf{x}_k),\\ \mathbf{k}_2 & = \mathbf{f}\!\left(t_k+\frac{h}{2},\mathbf{x}_k+\frac{h}{2}\mathbf{k}_1\right),\\ \mathbf{k}_3 & = \mathbf{f}\!\left(t_k+\frac{h}{2},\mathbf{x}_k+\frac{h}{2}\mathbf{k}_2\right),\\ \mathbf{k}_4 & = \mathbf{f}(t_k+h,\mathbf{x}_k+h\mathbf{k}_3),\\ \mathbf{x}_{k+1} & = \mathbf{x}_k + \frac{h}{6}\left(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4\right). \end{aligned} \]

flowchart TD
  I["Inputs: A(t), b(t), x0, step h, horizon N"] --> K["For k=0..N-1"]
  K --> F1["Compute f(tk, xk)"]
  F1 --> R["Choose method: Euler or RK4"]
  R -->|Euler| E["x(k+1) = x(k) + h*(A(tk)*x(k)+b(tk))"]
  R -->|RK4| G["k1,k2,k3,k4 then weighted update"]
  E --> U["Advance time: tk = tk + h"]
  G --> U
  U --> O["Store x(k+1) and continue"]
        

Practical note. For some matrices \( \mathbf{A} \), large \( h \) can cause numerical divergence even when the true solution behaves well. Step-size selection and stability will be revisited once we connect dynamics to eigen-structure more explicitly.

Standard scientific Python tools: numpy for linear algebra, scipy.integrate.solve_ivp for ODE integration. A commonly used control-focused library is python-control; we will use it more heavily once state-space models are formally introduced.

import numpy as np
from scipy.integrate import solve_ivp

# Example: time-varying linear system x_dot = A(t) x + b(t)
def A_of_t(t: float) -> np.ndarray:
    return np.array([[0.0, 1.0],
                     [-2.0 - 0.2*np.sin(t), -0.4]])

def b_of_t(t: float) -> np.ndarray:
    return np.array([0.0, 0.5*np.cos(t)])

def f(t: float, x: np.ndarray) -> np.ndarray:
    return A_of_t(t) @ x + b_of_t(t)

# Solve with SciPy
t0, tf = 0.0, 10.0
x0 = np.array([1.0, 0.0])
sol = solve_ivp(fun=f, t_span=(t0, tf), y0=x0, max_step=0.01, rtol=1e-8, atol=1e-10)

# From-scratch RK4 on the same system
def rk4_step(t: float, x: np.ndarray, h: float) -> np.ndarray:
    k1 = f(t, x)
    k2 = f(t + 0.5*h, x + 0.5*h*k1)
    k3 = f(t + 0.5*h, x + 0.5*h*k2)
    k4 = f(t + h,     x + h*k3)
    return x + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)

h = 0.01
N = int((tf - t0)/h)
ts = np.linspace(t0, tf, N+1)
xs = np.zeros((N+1, 2))
xs[0] = x0
t = t0
for k in range(N):
    xs[k+1] = rk4_step(t, xs[k], h)
    t += h

print("SciPy final state:", sol.y[:, -1])
print("RK4  final state:", xs[-1])
      

In modern control workflows, C++ is common for real-time simulation and embedded implementations. Use: Eigen for matrix algebra and either implement integrators directly or use Boost ODEInt.

#include <iostream>
#include <Eigen/Dense>
#include <cmath>

using Vec = Eigen::VectorXd;
using Mat = Eigen::MatrixXd;

Mat A_of_t(double t) {
    Mat A(2,2);
    A << 0.0, 1.0,
        -2.0 - 0.2*std::sin(t), -0.4;
    return A;
}

Vec b_of_t(double t) {
    Vec b(2);
    b << 0.0, 0.5*std::cos(t);
    return b;
}

Vec f(double t, const Vec& x) {
    return A_of_t(t)*x + b_of_t(t);
}

Vec rk4_step(double t, const Vec& x, double h) {
    Vec k1 = f(t, x);
    Vec k2 = f(t + 0.5*h, x + 0.5*h*k1);
    Vec k3 = f(t + 0.5*h, x + 0.5*h*k2);
    Vec k4 = f(t + h,     x + h*k3);
    return x + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
}

int main() {
    double t0 = 0.0, tf = 10.0, h = 0.01;
    int N = static_cast<int>((tf - t0)/h);

    Vec x(2);
    x << 1.0, 0.0;

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

    std::cout << "Final x = [" << x(0) << ", " << x(1) << "]\n";
    return 0;
}
      

In Java, a typical combination is EJML for linear algebra. For ODE integration, you can implement RK4 directly or use Apache Commons Math integrators (useful later in larger toolchains).

import org.ejml.simple.SimpleMatrix;

public class LinearODE_RK4 {
    static SimpleMatrix A_of_t(double t) {
        double a21 = -2.0 - 0.2*Math.sin(t);
        return new SimpleMatrix(new double[][] {
                {0.0, 1.0},
                {a21, -0.4}
        });
    }

    static SimpleMatrix b_of_t(double t) {
        return new SimpleMatrix(new double[][] {
                {0.0},
                {0.5*Math.cos(t)}
        });
    }

    static SimpleMatrix f(double t, SimpleMatrix x) {
        return A_of_t(t).mult(x).plus(b_of_t(t));
    }

    static SimpleMatrix rk4Step(double t, SimpleMatrix x, double h) {
        SimpleMatrix k1 = f(t, x);
        SimpleMatrix k2 = f(t + 0.5*h, x.plus(k1.scale(0.5*h)));
        SimpleMatrix k3 = f(t + 0.5*h, x.plus(k2.scale(0.5*h)));
        SimpleMatrix k4 = f(t + h,     x.plus(k3.scale(h)));
        return x.plus(k1.plus(k2.scale(2.0)).plus(k3.scale(2.0)).plus(k4).scale(h/6.0));
    }

    public static void main(String[] args) {
        double t0 = 0.0, tf = 10.0, h = 0.01;
        int N = (int)((tf - t0)/h);

        SimpleMatrix x = new SimpleMatrix(new double[][]{ {1.0},{0.0} });

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

        System.out.println("Final x = [" + x.get(0,0) + ", " + x.get(1,0) + "]");
    }
}
      

MATLAB provides robust variable-step integrators (e.g., ode45). In Simulink, the same dynamics can be built using integrator blocks and matrix gains. We keep the representation as “vector variables” here; later it becomes the standard internal form.

% Example: x_dot = A(t) x + b(t) with ode45
A_of_t = @(t) [0, 1; -2 - 0.2*sin(t), -0.4];
b_of_t = @(t) [0; 0.5*cos(t)];
f = @(t,x) A_of_t(t)*x + b_of_t(t);

tspan = [0, 10];
x0 = [1; 0];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,x] = ode45(f, tspan, x0, opts);

disp('Final x (ode45):');
disp(x(end,:).');
      

Programmatic Simulink construction (optional, code-driven). This creates a model that computes \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{b} \) using Integrator + Matrix Gain + Sum blocks for a constant \( \mathbf{A} \).

% Build a simple Simulink model for x_dot = A x + b (constant A, constant b)
A = [0 1; -2 -0.4];
b = [0; 0.5];

model = 'LinearVectorODE';
new_system(model); open_system(model);

add_block('simulink/Sources/Constant', [model '/b'], 'Value', 'b', 'Position', [30 80 80 110]);
add_block('simulink/Math Operations/Gain', [model '/A'], 'Gain', 'A', 'Multiplication', 'Matrix(K*u)', 'Position', [140 40 220 90]);
add_block('simulink/Math Operations/Sum', [model '/Sum'], 'Inputs', '++', 'Position', [270 55 295 85]);
add_block('simulink/Continuous/Integrator', [model '/Int'], 'Position', [340 40 370 90]);

% Initial condition
set_param([model '/Int'], 'InitialCondition', '[1;0]');

% Connections
add_line(model, 'Int/1', 'A/1');
add_line(model, 'A/1', 'Sum/1');
add_line(model, 'b/1', 'Sum/2');
add_line(model, 'Sum/1', 'Int/1');

set_param(model, 'StopTime', '10');
save_system(model);
      

In Simulink, the integrator output is the evolving vector \( \mathbf{x}(t) \), and the loop computes its derivative. Time-varying \( \mathbf{A}(t) \) or \( \mathbf{b}(t) \) can be implemented using function blocks, but the constant-matrix case is the clean baseline.

Mathematica supports direct vector ODE specification via NDSolve, and is also convenient for verifying numerical integrators.

(* Define A(t), b(t), and the vector ODE x'(t) = A(t) x(t) + b(t) *)
A[t_] := { {0, 1}, {-2 - 0.2 Sin[t], -0.4} };
b[t_] := {0, 0.5 Cos[t]};

eqns = {
  x1'[t] == (A[t].{x1[t], x2[t]})[[1]] + b[t][[1]],
  x2'[t] == (A[t].{x1[t], x2[t]})[[2]] + b[t][[2]],
  x1[0] == 1, x2[0] == 0
};

sol = NDSolve[eqns, {x1, x2}, {t, 0, 10}, AccuracyGoal -> 12, PrecisionGoal -> 12];
{x1[10], x2[10]} /. sol[[1]]
      

(* From-scratch RK4 for vector ODE *)
f[t_, x_] := A[t].x + b[t];

rk4Step[t_, x_, h_] := Module[{k1, k2, k3, k4},
  k1 = f[t, x];
  k2 = f[t + h/2, x + (h/2) k1];
  k3 = f[t + h/2, x + (h/2) k2];
  k4 = f[t + h,   x + h k3];
  x + (h/6) (k1 + 2 k2 + 2 k3 + k4)
];

h = 0.01; t0 = 0; tf = 10;
N = Round[(tf - t0)/h];
x = {1.0, 0.0}; t = t0;

Do[
  x = rk4Step[t, x, h];
  t = t + h;
, {k, 1, N}];

x
      

Problem 1 (Equivalence of scalar and vector forms): Let \( \mathbf{A}(t)=[a_{ij}(t)] \) and \( \mathbf{b}(t)=[b_i(t)] \). Prove that \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \) holds if and only if \( \dot{x}_i=\sum_{j=1}^n a_{ij}(t)x_j+b_i(t) \) for all \( i \).

Solution: Expand the matrix-vector product: \( \mathbf{A}(t)\mathbf{x}(t) \) has components \( (\mathbf{A}(t)\mathbf{x}(t))_i=\sum_{j=1}^n a_{ij}(t)x_j(t) \). Thus the vector equality is componentwise exactly the collection of scalar equalities, and conversely stacking scalar equalities yields the vector equality.

Problem 2 (Solution space is a vector space for homogeneous systems): Consider \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \) with continuous \( \mathbf{A}(t) \). Show that the set of solutions is closed under linear combinations.

Solution: If \( \mathbf{x}_1,\mathbf{x}_2 \) are solutions, then \( \dot{\mathbf{x}}_1=\mathbf{A}\mathbf{x}_1 \) and \( \dot{\mathbf{x}}_2=\mathbf{A}\mathbf{x}_2 \). For any scalars \( \alpha,\beta \):

\[ \frac{d}{dt}(\alpha\mathbf{x}_1+\beta\mathbf{x}_2) =\alpha\dot{\mathbf{x}}_1+\beta\dot{\mathbf{x}}_2 =\alpha\mathbf{A}\mathbf{x}_1+\beta\mathbf{A}\mathbf{x}_2 =\mathbf{A}(\alpha\mathbf{x}_1+\beta\mathbf{x}_2). \]

Hence \( \alpha\mathbf{x}_1+\beta\mathbf{x}_2 \) is also a solution.

Problem 3 (Affine structure of forced systems): Let \( \mathbf{x}_p \) be one solution of \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \). Prove that \( \mathbf{x} \) solves the forced system if and only if \( \mathbf{x}-\mathbf{x}_p \) solves the homogeneous system.

Solution: Define \( \mathbf{x}_h=\mathbf{x}-\mathbf{x}_p \). Then:

\[ \dot{\mathbf{x}}_h = \dot{\mathbf{x}}-\dot{\mathbf{x}}_p = (\mathbf{A}\mathbf{x}+\mathbf{b})-(\mathbf{A}\mathbf{x}_p+\mathbf{b}) = \mathbf{A}(\mathbf{x}-\mathbf{x}_p) = \mathbf{A}\mathbf{x}_h. \]

The reverse direction is identical by adding \( \mathbf{x}_p \) to any homogeneous solution.

Problem 4 (Effect of a constant similarity transformation): Let \( \mathbf{T} \) be invertible and \( \mathbf{z}=\mathbf{T}^{-1}\mathbf{x} \). Starting from \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \), derive the transformed system in \( \mathbf{z} \).

Solution: With constant \( \mathbf{T} \), \( \dot{\mathbf{x}}=\mathbf{T}\dot{\mathbf{z}} \). Substitute:

\[ \mathbf{T}\dot{\mathbf{z}}=\mathbf{A}(t)\mathbf{T}\mathbf{z}+\mathbf{b}(t) \;\Longrightarrow\; \dot{\mathbf{z}}=\mathbf{T}^{-1}\mathbf{A}(t)\mathbf{T}\mathbf{z}+\mathbf{T}^{-1}\mathbf{b}(t). \]

Problem 5 (One-step Euler update for constant coefficients): Suppose \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{b} \) where \( \mathbf{A},\mathbf{b} \) are constant. Derive the Euler recursion and show it is an affine map \( \mathbf{x}_{k+1}=\mathbf{M}\mathbf{x}_k+\mathbf{c} \).

Solution: Euler gives:

\[ \mathbf{x}_{k+1}=\mathbf{x}_k+h(\mathbf{A}\mathbf{x}_k+\mathbf{b}) =(\mathbf{I}+h\mathbf{A})\mathbf{x}_k + h\mathbf{b}. \]

Hence \( \mathbf{M}=\mathbf{I}+h\mathbf{A} \) and \( \mathbf{c}=h\mathbf{b} \).

We established the vector-matrix representation of coupled linear ODEs, developed the necessary calculus rules for vector and matrix functions, proved superposition properties for homogeneous systems and the affine structure of forced systems, and derived how coordinate changes transform linear dynamics via similarity. We also introduced Euler and RK4 time-stepping and implemented simulations in Python, C++, Java, MATLAB/Simulink, and Mathematica. The next lessons will build on this foundation to construct systematic solution operators (matrix exponential and fundamental matrix solutions).