Consider the linear (possibly time-varying) system in vector-matrix form:

\[ \dot{\mathbf{x}}(t) = \mathbf{A}(t)\mathbf{x}(t) + \mathbf{b}(t), \quad \mathbf{x}(t_0)=\mathbf{x}_0, \quad t \in [t_0, T] \]

where \( \mathbf{x}(t)\in\mathbb{R}^n \), \( \mathbf{A}(t)\in\mathbb{R}^{n\times n} \), and \( \mathbf{b}(t)\in\mathbb{R}^n \). The goal is to answer two questions:

• Existence: does at least one solution \( \mathbf{x}(\cdot) \) satisfy the IVP?
• Uniqueness: if a solution exists, is it the only one?

Throughout, we use a matrix norm \( \|\cdot\| \) compatible with the vector norm (submultiplicative): \( \|\mathbf{A}\mathbf{v}\| \le \|\mathbf{A}\|\|\mathbf{v}\| \).

If \( \mathbf{x}(\cdot) \) is differentiable and satisfies the differential equation, then integrating from \( t_0 \) to \( t \) yields the equivalent integral equation:

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

Conversely, if \( \mathbf{x}(\cdot) \) is continuous and satisfies the integral equation, then under mild regularity of \( \mathbf{A}(\cdot),\mathbf{b}(\cdot) \) it is differentiable and satisfies the ODE almost everywhere (and classically if the integrand is continuous). Thus, existence and uniqueness can be studied via a fixed-point problem on a space of continuous functions.

Define an operator \( \mathcal{T} \) on the space \( \mathcal{C}([t_0,t_0+h],\mathbb{R}^n) \) by:

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

A fixed point \( \mathbf{x}=\mathcal{T}\mathbf{x} \) is precisely a solution of the integral equation (hence of the IVP).

Assume \( \mathbf{A}(\cdot) \) and \( \mathbf{b}(\cdot) \) are continuous on \( [t_0,t_0+h] \). Then the quantities \( M := \max_{s\in [t_0,t_0+h]}\|\mathbf{A}(s)\| \) and \( B := \max_{s\in [t_0,t_0+h]}\|\mathbf{b}(s)\| \) are finite.

Contraction estimate. For any \( \mathbf{x}_1,\mathbf{x}_2 \in \mathcal{C} \),

\[ \|(\mathcal{T}\mathbf{x}_1)(t)-(\mathcal{T}\mathbf{x}_2)(t)\| \le \int_{t_0}^{t}\|\mathbf{A}(s)\|\cdot\|\mathbf{x}_1(s)-\mathbf{x}_2(s)\|\,ds. \]

Taking the supremum over \( t\in [t_0,t_0+h] \) gives, with the sup norm \( \|\mathbf{x}\|_{\infty}:=\sup_{t\in[t_0,t_0+h]}\|\mathbf{x}(t)\| \),

\[ \|\mathcal{T}\mathbf{x}_1-\mathcal{T}\mathbf{x}_2\|_{\infty} \le Mh\;\|\mathbf{x}_1-\mathbf{x}_2\|_{\infty}. \]

Therefore, if \( Mh < 1 \), then \( \mathcal{T} \) is a contraction on \( \mathcal{C}([t_0,t_0+h],\mathbb{R}^n) \). By the Banach fixed-point theorem, there exists a unique fixed point, which is the unique solution on \( [t_0,t_0+h] \).

Picard iteration (constructive existence). Start with \( \mathbf{x}^{(0)}(t)=\mathbf{x}_0 \) and iterate

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

Under \( Mh < 1 \), \( \mathbf{x}^{(k)} \) converges uniformly to the unique solution.

flowchart TD
  A["Given A(t), b(t), x0, interval [t0,T]"] --> B["Convert to integral equation x = T(x)"]
  B --> C["Choose step h so that M*h < 1 (M = max ||A(t)||)"]
  C --> D["Run Picard iteration: x(k+1)=T(x(k))"]
  D --> E["Converges to unique fixed point on [t0,t0+h]"]
  E --> F["Extend interval by repeating steps to reach T"]
        

The contraction argument is local because it requires \( Mh < 1 \). However, if \( \mathbf{A}(\cdot) \) and \( \mathbf{b}(\cdot) \) are continuous on the entire finite interval \( [t_0,T] \), then \( \max_{t\in[t_0,T]}\|\mathbf{A}(t)\| \) is finite. We can partition \( [t_0,T] \) into subintervals of length \( h \) small enough to satisfy the contraction condition on each subinterval and “march” the unique solution forward in time.

This yields:

Theorem (Existence and Uniqueness on a finite interval).

If \( \mathbf{A}(\cdot) \) and \( \mathbf{b}(\cdot) \) are continuous on \( [t_0,T] \), then for every initial condition \( \mathbf{x}(t_0)=\mathbf{x}_0 \) there exists a unique solution \( \mathbf{x}(\cdot) \) on \( [t_0,T] \).

Because the system is linear (hence globally Lipschitz in \( \mathbf{x} \) on bounded time intervals), no finite-time blow-up occurs on any finite interval when \( \mathbf{A}(t) \) is bounded.

Let \( \mathbf{x}_1(\cdot) \) and \( \mathbf{x}_2(\cdot) \) be two solutions of the same IVP (same \( \mathbf{A}(t),\mathbf{b}(t),\mathbf{x}_0 \)). Define the difference \( \mathbf{e}(t):=\mathbf{x}_1(t)-\mathbf{x}_2(t) \). Then:

\[ \dot{\mathbf{e}}(t) = \mathbf{A}(t)\mathbf{e}(t), \quad \mathbf{e}(t_0)=\mathbf{0}. \]

Integrating,

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

Taking norms and using submultiplicativity:

\[ \|\mathbf{e}(t)\| \le \int_{t_0}^{t}\|\mathbf{A}(s)\|\cdot \|\mathbf{e}(s)\|\,ds. \]

A standard Grönwall inequality argument implies \( \|\mathbf{e}(t)\|=0 \) for all \( t\in[t_0,T] \), hence \( \mathbf{x}_1(t)=\mathbf{x}_2(t) \) for all \( t \). This proves uniqueness.

Continuous dependence on initial data. If \( \mathbf{x}_1,\mathbf{x}_2 \) solve the same dynamics but with different initial conditions \( \mathbf{x}_1(t_0)=\mathbf{x}_{0,1} \), \( \mathbf{x}_2(t_0)=\mathbf{x}_{0,2} \), then

\[ \|\mathbf{x}_1(t)-\mathbf{x}_2(t)\| \le \|\mathbf{x}_{0,1}-\mathbf{x}_{0,2}\|\; \exp\!\left(\int_{t_0}^{t}\|\mathbf{A}(s)\|\,ds\right). \]

This estimate is foundational: it shows trajectories are not only unique, but also stable (in a mathematical sense) with respect to small perturbations in initial conditions when \( \int_{t_0}^{t}\|\mathbf{A}(s)\|ds \) is moderate.

From Lesson 4, the homogeneous system \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \) admits a fundamental matrix \( \boldsymbol{\Phi}(t,t_0) \) satisfying:

\[ \frac{d}{dt}\boldsymbol{\Phi}(t,t_0)=\mathbf{A}(t)\boldsymbol{\Phi}(t,t_0), \quad \boldsymbol{\Phi}(t_0,t_0)=\mathbf{I}. \]

Existence and uniqueness for the matrix IVP above follows by applying the same theorem componentwise (or columnwise): each column solves a linear IVP of the form \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \), hence is uniquely defined.

For the nonhomogeneous system, the unique solution is given by the variation-of-constants formula:

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

In the LTI case \( \mathbf{A}(t)=\mathbf{A} \) constant, this reduces to the matrix exponential representation:

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

Numerical solvers (Euler, Runge–Kutta, etc.) approximate a trajectory that is assumed to be well-defined. Existence ensures there is a trajectory to approximate; uniqueness ensures that different numerical methods (with refinement) converge toward the same physical trajectory rather than multiple incompatible ones.

flowchart TD
  A["Model: xdot = A(t)x + b(t), x(t0)=x0"] --> B["Check conditions: A(t), b(t) continuous on [t0,T]"]
  B --> C["Well-posed IVP: unique trajectory exists"]
  C --> D["Pick numerical method: Euler / RK / expm-based"]
  D --> E["Refine step size and compare solvers"]
  E --> F["Agreement indicates convergence to the unique solution"]
        

We simulate a time-varying linear system and empirically verify uniqueness by comparing (i) a high-order adaptive integrator against (ii) a from-scratch Euler method as the step size decreases.

import numpy as np

# Time-varying linear system: xdot = A(t)x + b(t)
def A(t):
    # bounded, continuous A(t) on any finite interval
    return np.array([[0.0, 1.0],
                     [-2.0 - 0.5*np.sin(t), -0.4]])

def b(t):
    return np.array([0.0, 1.0*np.cos(2.0*t)])

def f(t, x):
    return A(t) @ x + b(t)

t0, T = 0.0, 8.0
x0 = np.array([1.0, 0.0])

# 1) Reference solution with SciPy's solve_ivp (adaptive RK method)
from scipy.integrate import solve_ivp
sol = solve_ivp(lambda t, x: f(t, x), (t0, T), x0, method="RK45", rtol=1e-9, atol=1e-12)
t_ref = sol.t
x_ref = sol.y.T

# 2) From-scratch explicit Euler method with uniform step
def euler_uniform(f, t0, T, x0, N):
    h = (T - t0) / N
    t = t0
    x = x0.astype(float).copy()
    traj_t = [t]
    traj_x = [x.copy()]
    for k in range(N):
        x = x + h * f(t, x)
        t = t + h
        traj_t.append(t)
        traj_x.append(x.copy())
    return np.array(traj_t), np.array(traj_x)

# Compare Euler with decreasing step sizes against reference (interpolate reference)
from scipy.interpolate import interp1d
x_ref_interp = interp1d(t_ref, x_ref, axis=0, kind="cubic", fill_value="extrapolate")

for N in [200, 400, 800, 1600]:
    t_eu, x_eu = euler_uniform(f, t0, T, x0, N)
    err = np.max(np.linalg.norm(x_eu - x_ref_interp(t_eu), axis=1))
    print(f"N={N:4d}, h={(T-t0)/N:.5f}, max error ~ {err:.3e}")

# Uniqueness sanity check: same IVP should yield same trajectory (within numerical tolerances)
sol2 = solve_ivp(lambda t, x: f(t, x), (t0, T), x0, method="DOP853", rtol=1e-10, atol=1e-13)
x2_interp = interp1d(sol2.t, sol2.y.T, axis=0, kind="cubic", fill_value="extrapolate")
grid = np.linspace(t0, T, 401)
diff = np.max(np.linalg.norm(x_ref_interp(grid) - x2_interp(grid), axis=1))
print("Max difference between two high-accuracy solvers:", diff)
      

Practical interpretation: if \( \mathbf{A}(t) \) and \( \mathbf{b}(t) \) are continuous on \( [t_0,T] \), the IVP is well-posed and the two high-accuracy solvers should agree up to numerical tolerance; Euler converges as \( h \) decreases because it is approximating a unique trajectory.

Relevant Python libraries commonly used later in modern control workflows include: numpy, scipy (ODE solvers and linear algebra), and the control toolbox (state-space objects and responses). At this chapter, we mainly rely on numerical integration and matrix operations.

The following example integrates \( \dot{x}=A(t)x+b(t) \) using Boost.Odeint and Eigen for vector/matrix operations.

#include <iostream>
#include <cmath>
#include <Eigen/Dense>
#include <boost/numeric/odeint.hpp>

using Vector = Eigen::Vector2d;
using Matrix = Eigen::Matrix2d;

struct LinSystem {
  void operator()(const Vector &x, Vector &dxdt, const double t) const {
    Matrix A;
    A << 0.0, 1.0,
         -2.0 - 0.5*std::sin(t), -0.4;
    Vector b;
    b << 0.0, std::cos(2.0*t);
    dxdt = A * x + b;
  }
};

int main() {
  using namespace boost::numeric::odeint;

  Vector x;
  x << 1.0, 0.0;

  double t0 = 0.0, T = 8.0, dt = 1e-3;

  auto obs = [](const Vector &x, double t) {
    if (std::fmod(t, 1.0) < 1e-12) {
      std::cout << "t=" << t << "  x=[" << x(0) << ", " << x(1) << "]\n";
    }
  };

  runge_kutta_dopri5<Vector, double, Vector, double, vector_space_algebra> stepper;
  integrate_const(stepper, LinSystem{}, x, t0, T, dt, obs);

  return 0;
}
      

Notes: Continuity and boundedness of \( A(t) \) on finite intervals is a theoretical guarantee for a unique solution, while numerical integrators approximate that unique trajectory. Eigen provides efficient linear algebra; Odeint provides robust integrators suitable for control simulation pipelines.

Java users often rely on Apache Commons Math for ODE integration. The example below implements \( \dot{x}=A(t)x+b(t) \) via an ODE interface and uses a Dormand–Prince integrator.

import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.ode.sampling.StepHandler;
import org.apache.commons.math3.ode.sampling.StepInterpolator;

public class LinearTimeVaryingODE {

  static class LinSys implements FirstOrderDifferentialEquations {
    @Override
    public int getDimension() { return 2; }

    @Override
    public void computeDerivatives(double t, double[] x, double[] dxdt) {
      double a00 = 0.0, a01 = 1.0;
      double a10 = -2.0 - 0.5*Math.sin(t);
      double a11 = -0.4;

      double b0 = 0.0;
      double b1 = Math.cos(2.0*t);

      dxdt[0] = a00*x[0] + a01*x[1] + b0;
      dxdt[1] = a10*x[0] + a11*x[1] + b1;
    }
  }

  public static void main(String[] args) {
    double t0 = 0.0, T = 8.0;
    double[] x = new double[] { 1.0, 0.0 };

    DormandPrince853Integrator integrator =
        new DormandPrince853Integrator(1e-6, 1e-2, 1e-10, 1e-12);

    integrator.addStepHandler(new StepHandler() {
      @Override
      public void init(double t0, double[] x0, double t) { }

      @Override
      public void handleStep(StepInterpolator interpolator, boolean isLast) {
        double t = interpolator.getCurrentTime();
        if (Math.abs(t - Math.rint(t)) < 1e-3) { // near integer times
          double[] state = interpolator.getInterpolatedState();
          System.out.printf("t=%.3f  x=[%.6f, %.6f]%n", t, state[0], state[1]);
        }
      }
    });

    integrator.integrate(new LinSys(), t0, x, T, x);
    System.out.printf("Final state x(T)=[%.8f, %.8f]%n", x[0], x[1]);
  }
}
      

In later modern-control implementations (observers, estimators, discrete-time controllers), you will often combine: (i) linear algebra libraries (EJML, Apache Commons Math linear) and (ii) numerical ODE solvers. The theoretical uniqueness result ensures the model you simulate has a single consistent trajectory.

MATLAB (ODE solver + fundamental matrix idea).

% Time-varying linear system: xdot = A(t)x + b(t)
A = @(t) [0, 1;
          -2 - 0.5*sin(t), -0.4];
b = @(t) [0; cos(2*t)];
f = @(t,x) A(t)*x + b(t);

t0 = 0; T = 8;
x0 = [1; 0];

opts = odeset('RelTol',1e-9,'AbsTol',1e-12);
[t,x] = ode45(f, [t0 T], x0, opts);

disp(['x(T) = [', num2str(x(end,1)), ', ', num2str(x(end,2)), ']']);

% LTI special case check (when A is constant):
A0 = [0 1; -2 -0.4];
x_lti = expm(A0*(T-t0))*x0;  % homogeneous solution for b(t)=0
      

Simulink (state-space block).

• Use a State-Space block for constant \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \).
• For time-varying \( \mathbf{A}(t) \) or \( \mathbf{b}(t) \), use an Integrator block with a MATLAB Function block computing \( \dot{x} \).
• Set solver to variable-step (e.g., ode45 or ode15s if stiffness appears).

Wolfram Mathematica (symbolic/numeric integration).

(* Time-varying linear system: x'(t) = A(t) x(t) + b(t) *)
A[t_] := { {0, 1}, {-2 - 0.5 Sin[t], -0.4} };
b[t_] := {0, Cos[2 t]};

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

sol = NDSolve[eqns, {x1, x2}, {t, 0, 8}, WorkingPrecision -> 40][[1]];
{x1[8], x2[8]} /. sol

(* LTI case (constant A) via MatrixExp *)
A0 = { {0, 1}, {-2, -0.4} };
x0 = {1, 0};
MatrixExp[A0*8].x0
      

Mathematica offers high-precision ODE integration and closed-form manipulations; MATLAB/Simulink offers standard engineering simulation workflows. In both, the correctness of the simulation fundamentally relies on well-posedness (existence + uniqueness) of the IVP.

Problem 1 (Lipschitz constant for linear dynamics): Let \( f(t,\mathbf{x})=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \) with \( \mathbf{A}(t) \) continuous on \( [t_0,T] \). Show that \( f \) is Lipschitz in \( \mathbf{x} \) uniformly in \( t \) on \( [t_0,T] \), and identify a valid Lipschitz constant.

Solution: For any \( \mathbf{x}_1,\mathbf{x}_2 \),

\[ \|f(t,\mathbf{x}_1)-f(t,\mathbf{x}_2)\| = \|\mathbf{A}(t)(\mathbf{x}_1-\mathbf{x}_2)\| \le \|\mathbf{A}(t)\|\;\|\mathbf{x}_1-\mathbf{x}_2\|. \]

Let \( L := \max_{t\in[t_0,T]}\|\mathbf{A}(t)\| \) (finite by continuity on a compact interval). Then \( \|f(t,\mathbf{x}_1)-f(t,\mathbf{x}_2)\| \le L\|\mathbf{x}_1-\mathbf{x}_2\| \) for all \( t\in[t_0,T] \). Thus \( L \) is a Lipschitz constant.

Problem 2 (Uniqueness by Grönwall): Suppose \( \mathbf{x}_1,\mathbf{x}_2 \) both solve the same IVP \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \), \( \mathbf{x}(t_0)=\mathbf{x}_0 \). Prove \( \mathbf{x}_1(t)=\mathbf{x}_2(t) \) for all \( t\in[t_0,T] \).

Solution: Let \( \mathbf{e}(t)=\mathbf{x}_1(t)-\mathbf{x}_2(t) \). Then

\[ \dot{\mathbf{e}}(t)=\mathbf{A}(t)\mathbf{e}(t), \quad \mathbf{e}(t_0)=\mathbf{0}. \]

Integrating gives \( \mathbf{e}(t)=\int_{t_0}^{t}\mathbf{A}(s)\mathbf{e}(s)\,ds \), hence

\[ \|\mathbf{e}(t)\| \le \int_{t_0}^{t}\|\mathbf{A}(s)\|\;\|\mathbf{e}(s)\|\,ds. \]

By Grönwall, the only continuous function satisfying this inequality with \( \mathbf{e}(t_0)=0 \) is \( \mathbf{e}(t)\equiv 0 \). Therefore \( \mathbf{x}_1(t)=\mathbf{x}_2(t) \) for all \( t\in[t_0,T] \).

Problem 3 (Bound on solutions): Assume \( \mathbf{A}(t),\mathbf{b}(t) \) are continuous on \( [t_0,T] \). Show that any solution satisfies an a priori bound of the form \( \|\mathbf{x}(t)\| \le \alpha(t) \) for an explicit function \( \alpha(t) \) depending on \( \|\mathbf{x}_0\| \), \( \max\|\mathbf{A}(t)\| \), and \( \max\|\mathbf{b}(t)\| \).

Solution: From the integral equation,

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

Let \( L:=\max_{s\in[t_0,T]}\|\mathbf{A}(s)\| \) and \( B:=\max_{s\in[t_0,T]}\|\mathbf{b}(s)\| \). Then

\[ \|\mathbf{x}(t)\| \le \|\mathbf{x}_0\| + B(t-t_0) + \int_{t_0}^{t}L\|\mathbf{x}(s)\|\,ds. \]

Applying Grönwall to the function \( u(t)=\|\mathbf{x}(t)\| \) yields

\[ \|\mathbf{x}(t)\| \le \left(\|\mathbf{x}_0\| + B(t-t_0)\right)\exp\!\left(L(t-t_0)\right), \quad t\in[t_0,T]. \]

This bound is useful for showing that solutions cannot blow up on finite time intervals when \( L \) is finite.

Problem 4 (LTI uniqueness implies identical state transition): Consider the LTI homogeneous system \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \). Let \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \). Prove that any other solution with the same initial condition must coincide with this \( \mathbf{x}(t) \).

Solution: Let \( \mathbf{y}(t) \) be any solution with \( \mathbf{y}(t_0)=\mathbf{x}_0 \), and define \( \mathbf{e}(t)=\mathbf{y}(t)-e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \). Then \( \dot{\mathbf{e}}=\mathbf{A}\mathbf{e} \) and \( \mathbf{e}(t_0)=0 \). By the uniqueness result for linear systems (Problem 2 with constant \( \mathbf{A}(t)=\mathbf{A} \), \( \mathbf{b}(t)=0 \)), \( \mathbf{e}(t)\equiv 0 \), hence \( \mathbf{y}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \).

Problem 5 (Stepwise extension argument): Suppose \( \mathbf{A}(t) \) is continuous on \( [t_0,T] \). Explain how to choose a uniform step length \( h \) so that the contraction condition holds on each subinterval and show that the unique local solutions can be concatenated to obtain a unique global solution on \( [t_0,T] \).

Solution: Since \( \mathbf{A}(t) \) is continuous on a compact interval, \( L=\max_{t\in[t_0,T]}\|\mathbf{A}(t)\| < \infty \). Choose \( h \) such that \( Lh < 1 \) (e.g., \( h = \min\{(1/2L),\, T-t_0\} \) if \( L>0 \), and any \( h \) if \( L=0 \)). Partition \( [t_0,T] \) into subintervals of length at most \( h \). On each subinterval, the Picard operator is a contraction, giving a unique solution there. The endpoint of one subinterval provides the initial condition for the next; uniqueness on overlaps ensures the pieces match at boundaries, producing a single continuous (indeed differentiable) solution on the entire interval.

We proved that the linear IVP \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x}+\mathbf{b}(t) \), \( \mathbf{x}(t_0)=\mathbf{x}_0 \), is well-posed on any finite interval when \( \mathbf{A}(t) \) and \( \mathbf{b}(t) \) are continuous. The key steps were: (i) rewrite as an integral equation, (ii) use contraction mapping (Picard iteration) for local existence and uniqueness, (iii) extend stepwise to global solutions on \( [t_0,T] \), and (iv) establish uniqueness and continuous dependence via Grönwall-type estimates. These results provide the rigorous foundation for using fundamental matrices and matrix exponentials as solution operators in later chapters.