[proofplan] We construct the principal fundamental matrix column by column, taking the $j$-th column to be the solution with initial value the $j$-th standard basis vector. The local input used below is the already staged maximal-solution theorem [citetheorem:8375]: for a continuous vector field on $I\times E$ that is locally Lipschitz in the state variable locally uniformly in time, and for any initial condition in $I\times E$, it gives a unique maximal solution on an interval open in the relative topology of $I$, with the maximal restriction property for every other solution of the same [initial value problem](/page/Initial%20Value%20Problem). The only point requiring care is that the linear equation is solved on the whole interval $I$, not merely on this maximal interval; a compact-interval boundedness estimate prevents finite-time blow-up and forces the maximal interval to be all of $I$. Once the columns are constructed, the matrix equation and normalization are immediate. Invertibility at every time follows from uniqueness: a non-trivial linear dependence at time $t_*$ would give a solution that vanishes at $t_*$ and hence vanishes identically, contradicting the normalization at $t_0$. Finally, uniqueness of the matrix and the representation $x(t)=\Phi(t)x_0$ both follow by applying uniqueness to vector-valued initial value problems. [/proofplan]

[step:Obtain global solutions for all initial vectors in $\mathbb{R}^n$]Let $b \in \mathbb{R}^n$. Throughout this proof, $\mathcal{L}^1$ denotes one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}$, and $\|B\|_{\mathrm{op}}$ denotes the operator norm of a matrix $B \in \mathbb{R}^{n \times n}$ acting on $\mathbb{R}^n$ with the Euclidean norm. Define the time-dependent vector field $f: I \times \mathbb{R}^n \to \mathbb{R}^n$ by $(t,x) \mapsto A(t)x$. Since $A$ is continuous, $f$ is continuous. Moreover, on each compact subinterval $K \subset I$, the quantity \begin{align*} M_K := \sup_{t \in K} \|A(t)\|_{\mathrm{op}} \end{align*} is finite, and for all $t \in K$ and all $x,y \in \mathbb{R}^n$ one has \begin{align*} |f(t,x)-f(t,y)| \leq M_K |x-y|. \end{align*} Thus $f$ is locally Lipschitz in the state variable locally uniformly in time, with the compact interval $K$ supplying a uniform Lipschitz constant on each local time neighbourhood. The theorem [citetheorem:8375] is used here in its formalized maximal-solution form: its hypotheses are exactly the continuity of $f$ on $I\times E$, openness of $E$, and local uniform state-Lipschitz condition, and its conclusion supplies a unique maximal solution on a relatively open interval of $I$ together with the maximal restriction property. Applying it with $E=\mathbb{R}^n$ gives a unique maximal solution \begin{align*} u_b: J_b \to \mathbb{R}^n \end{align*} of $\dot{u}_b(t)=A(t)u_b(t)$ with $u_b(t_0)=b$, where $J_b \subset I$ is an interval containing $t_0$ and open in the relative topology of $I$. The relative-topology formulation includes the case in which the initial time is an endpoint of $I$. We prove that $J_b=I$. Let $K \subset I$ be a compact interval with $t_0 \in K$. For $t \in J_b \cap K$, the integral form of the differential equation gives \begin{align*} u_b(t)=b+\int_{t_0}^{t} A(s)u_b(s)\, d\mathcal{L}^1(s), \end{align*} where the integral over an oriented interval is understood with the usual sign convention. Taking Euclidean norms and using the operator norm gives \begin{align*} |u_b(t)| \leq |b|+\int_{\min\{t,t_0\}}^{\max\{t,t_0\}} M_K |u_b(s)|\, d\mathcal{L}^1(s). \end{align*} The integral Gronwall estimate in the following form applies on the compact interval with endpoints $t_0$ and $t$: if a non-negative [continuous function](/page/Continuous%20Function) $g$ satisfies \begin{align*} g(r) \leq a+\int_{\min\{r,t_0\}}^{\max\{r,t_0\}} M_K g(s)\, d\mathcal{L}^1(s) \end{align*} for $r$ between $t_0$ and $t$, then $g(t) \leq a e^{M_K |t-t_0|}$. Applying this with $g(r)=|u_b(r)|$ and $a=|b|$ yields \begin{align*} |u_b(t)| \leq |b| e^{M_K |t-t_0|}. \end{align*} Hence $u_b$ is bounded on $J_b \cap K$. Suppose for contradiction that $J_b \neq I$. Since $J_b$ is relatively open in the interval $I$, there is an endpoint $\tau \in I$ of $J_b$ that belongs to the closure of $J_b$ but not to $J_b$. Choose a compact interval $K \subset I$ containing $t_0$ and $\tau$. The preceding bound gives $|u_b(t)| \leq C_K$ for all $t \in J_b \cap K$, where $d_K := \max K - \min K$ denotes the length of the compact interval $K$ and \begin{align*} C_K := |b| e^{M_K d_K}. \end{align*} For $r,s \in J_b \cap K$, the integral equation gives \begin{align*} |u_b(r)-u_b(s)| \leq \int_{\min\{r,s\}}^{\max\{r,s\}} M_K C_K\, d\mathcal{L}^1(\theta)=M_K C_K |r-s|. \end{align*} Thus $u_b(t)$ has a finite limit $b_\tau \in \mathbb{R}^n$ as $t \to \tau$ with $t \in J_b$. Apply [citetheorem:8375] again with initial data $(\tau,b_\tau)$. Since the theorem is stated for intervals and neighbourhoods relative to $I$, this is valid even if $\tau$ is an endpoint of $I$. Its maximal solution contains a relatively open neighbourhood of $\tau$ in $I$, and restricting that maximal solution to such a neighbourhood gives a local solution \begin{align*} w: N \to \mathbb{R}^n \end{align*} through $(\tau,b_\tau)$ on a relative neighbourhood $N \subset I$ of $\tau$. We verify compatibility on the side where $N$ overlaps $J_b$. Let $r \in N \cap J_b$. Since $u_b(t) \to b_\tau$ as $t \to \tau$ through $J_b$, and since $A(s)u_b(s)$ is bounded on $J_b \cap K$, absolute continuity of the [Lebesgue integral](/page/Lebesgue%20Integral) on the compact interval $K$ lets the integral identity for $u_b$ between $r$ and $t \in J_b$ pass to the limit $t \to \tau$ and gives \begin{align*} u_b(r)=b_\tau+\int_\tau^r A(s)u_b(s)\, d\mathcal{L}^1(s), \end{align*} with the oriented-interval convention. Thus the continuous extension of $u_b$ to $\tau$ solves the same initial value problem as $w$ on the one-sided interval $N \cap J_b \cup \{\tau\}$. The local uniqueness consequence contained in [citetheorem:8375], applied at the initial condition $(\tau,b_\tau)$ and then restricted to the common side of $\tau$, gives $w(t)=u_b(t)$ for $t \in N \cap J_b$. Hence the piecewise map obtained by attaching $w$ to $u_b$ at $\tau$ is a solution on a strictly larger interval in $I$. This contradicts maximality. Therefore no such endpoint exists, and $J_b=I$.[/step]

[guided]We first explain why the solution exists on all of $I$. Here $\mathcal{L}^1$ denotes one-dimensional Lebesgue measure on $\mathbb{R}$, and $\|B\|_{\mathrm{op}}$ denotes the operator norm of a matrix $B \in \mathbb{R}^{n \times n}$ acting on $\mathbb{R}^n$ with the Euclidean norm. The cited theorem, Existence and Uniqueness of the Maximal Solution [citetheorem:8375], applies to the map \begin{align*} f: I \times \mathbb{R}^n \to \mathbb{R}^n, \qquad (t,x) \mapsto A(t)x. \end{align*} Indeed, $f$ is continuous because $A$ is continuous and matrix-vector multiplication is continuous. To check the local Lipschitz condition in the state variable, fix a compact interval $K \subset I$ and define \begin{align*} M_K := \sup_{t \in K} \|A(t)\|_{\mathrm{op}}. \end{align*} The supremum is finite because $A$ is continuous on the compact set $K$. For $t \in K$ and $x,y \in \mathbb{R}^n$, \begin{align*} |f(t,x)-f(t,y)|=|A(t)(x-y)| \leq \|A(t)\|_{\mathrm{op}} |x-y| \leq M_K |x-y|. \end{align*} Thus the hypotheses of this maximal-solution theorem are satisfied, so for each $b \in \mathbb{R}^n$ there is a unique maximal solution \begin{align*} u_b: J_b \to \mathbb{R}^n \end{align*} with $u_b(t_0)=b$. The remaining issue is global existence on $I$. A local existence theorem by itself only gives a maximal interval $J_b$. We show that the solution cannot break down at a finite point of $I$. Let $K \subset I$ be a compact interval containing $t_0$. For $t \in J_b \cap K$, the differential equation is equivalent to the integral identity \begin{align*} u_b(t)=b+\int_{t_0}^{t} A(s)u_b(s)\, d\mathcal{L}^1(s). \end{align*} Taking Euclidean norms and estimating the integrand by the operator norm gives \begin{align*} |u_b(t)| \leq |b|+\int_{\min\{t,t_0\}}^{\max\{t,t_0\}} M_K |u_b(s)|\, d\mathcal{L}^1(s). \end{align*} We use the integral Gronwall estimate in its orientation-free compact-interval form: if a non-negative continuous function $g$ satisfies \begin{align*} g(r) \leq a+\int_{\min\{r,t_0\}}^{\max\{r,t_0\}} M_K g(s)\, d\mathcal{L}^1(s) \end{align*} for $r$ between $t_0$ and $t$, then $g(t) \leq a e^{M_K |t-t_0|}$. Applying this to $g(r)=|u_b(r)|$ with $a=|b|$ gives \begin{align*} |u_b(t)| \leq |b| e^{M_K |t-t_0|}. \end{align*} This is the key bound: on every compact subinterval of $I$, the solution remains bounded. Now suppose $J_b$ were not all of $I$. Since $J_b$ is relatively open in the interval $I$, it would have an endpoint $\tau \in I$ lying in the closure of $J_b$ but not in $J_b$. Choose a compact interval $K \subset I$ containing both $t_0$ and $\tau$. The bound above gives \begin{align*} |u_b(t)| \leq C_K \end{align*} for all $t \in J_b \cap K$, where $d_K := \max K - \min K$ denotes the length of the compact interval $K$ and \begin{align*} C_K := |b| e^{M_K d_K}. \end{align*} For $r,s \in J_b \cap K$, using the integral equation between $r$ and $s$ gives \begin{align*} |u_b(r)-u_b(s)| \leq \int_{\min\{r,s\}}^{\max\{r,s\}} |A(\theta)u_b(\theta)|\, d\mathcal{L}^1(\theta) \leq M_K C_K |r-s|. \end{align*} Thus $u_b$ is Cauchy as $t \to \tau$ through points of $J_b$, and therefore has a limit $b_\tau \in \mathbb{R}^n$. Finally apply the maximal-solution theorem [citetheorem:8375] with initial data $(\tau,b_\tau)$. The theorem uses intervals and neighbourhoods relative to $I$, so it applies at endpoint times as well as interior times. Its maximal solution restricts to a local solution \begin{align*} w: N \to \mathbb{R}^n \end{align*} starting at time $\tau$ with value $b_\tau$ on a relative neighbourhood $N \subset I$ of $\tau$. We must first check that the old solution really has initial value $b_\tau$ at $\tau$ in the integral sense. Fix $r \in N \cap J_b$. For $t \in J_b$ on the same side of $\tau$ as $r$, the integral equation gives \begin{align*} u_b(r)=u_b(t)+\int_t^r A(s)u_b(s)\, d\mathcal{L}^1(s). \end{align*} Letting $t \to \tau$ through points of $J_b$, the first term tends to $b_\tau$. The integrands $A(s)u_b(s)$ are bounded by $M_KC_K$ on $J_b \cap K$, so absolute continuity of the Lebesgue integral on $K$ implies that the integral over the shrinking interval from $t$ to $\tau$ tends to $0$, and we obtain \begin{align*} u_b(r)=b_\tau+\int_\tau^r A(s)u_b(s)\, d\mathcal{L}^1(s). \end{align*} Thus the continuous extension of $u_b$ to $\tau$ solves the same initial value problem through $(\tau,b_\tau)$ on the one-sided overlap. The uniqueness clause of [citetheorem:8375] therefore forces $w$ to agree with this extension on $N \cap J_b$. Hence $w$ can be joined to $u_b$ at $\tau$ to form a solution on a strictly larger interval in $I$. This contradicts maximality of $u_b$. Hence no endpoint obstruction exists, and the maximal domain is $J_b=I$.[/guided]

[step:Build the matrix from the standard basis solutions] For each $j \in \{1,\dots,n\}$, let $e_j \in \mathbb{R}^n$ denote the $j$-th standard basis vector. From the previous step, let \begin{align*} u_j: I \to \mathbb{R}^n \end{align*} be the unique solution of \begin{align*} \dot{u}_j(t)=A(t)u_j(t), \qquad u_j(t_0)=e_j. \end{align*} Define \begin{align*} \Phi: I \to \mathbb{R}^{n \times n} \end{align*} by declaring its $j$-th column to be $u_j(t)$ for each $j \in \{1,\dots,n\}$ and each $t \in I$. Because each $u_j$ is continuous on $I$ and differentiable on the interior of $I$, the matrix-valued map $\Phi$ is continuous on $I$ and differentiable on the interior of $I$. For every interior point $t$ of $I$, the $j$-th column of $\Phi'(t)$ is $\dot{u}_j(t)$, while the $j$-th column of $A(t)\Phi(t)$ is $A(t)u_j(t)$. Since $\dot{u}_j(t)=A(t)u_j(t)$ for every $j$, we obtain \begin{align*} \Phi'(t)=A(t)\Phi(t). \end{align*} At $t=t_0$, the $j$-th column of $\Phi(t_0)$ is $e_j$, hence \begin{align*} \Phi(t_0)=I_n. \end{align*} It remains to prove that $\Phi(t)$ is invertible for every $t \in I$. Fix $t_* \in I$, and suppose $c \in \mathbb{R}^n$ satisfies $\Phi(t_*)c=0$. Define \begin{align*} z: I \to \mathbb{R}^n, \qquad t \mapsto \Phi(t)c. \end{align*} Then $z$ is continuous on $I$, differentiable on the interior of $I$, and satisfies \begin{align*} \dot{z}(t)=A(t)z(t) \end{align*} for every interior point $t$ of $I$. Also $z(t_*)=0$. The zero map \begin{align*} 0_I: I \to \mathbb{R}^n, \qquad t \mapsto 0, \end{align*} solves the same initial value problem at time $t_*$. The hypotheses of the maximal-solution theorem [citetheorem:8375] have already been verified for the vector field $(t,x) \mapsto A(t)x$. Since that theorem permits initial times anywhere in $I$ with domains open in the relative topology of $I$, it applies to the initial time $t_*$. Its uniqueness clause gives $z=0_I$ on $I$. Evaluating at $t_0$ gives \begin{align*} c=I_nc=\Phi(t_0)c=z(t_0)=0. \end{align*} Therefore $\ker \Phi(t_*)=\{0\}$. Since $\Phi(t_*)$ is an $n \times n$ real matrix, injectivity implies invertibility. Because $t_* \in I$ was arbitrary, $\Phi(t)$ is invertible for every $t \in I$. Thus $\Phi$ is a principal fundamental matrix based at $t_0$. [/step]

[step:Prove uniqueness of the principal fundamental matrix column by column] Let \begin{align*} \Psi: I \to \mathbb{R}^{n \times n} \end{align*} be another continuous matrix-valued map, differentiable on the interior of $I$, satisfying \begin{align*} \Psi'(t)=A(t)\Psi(t) \end{align*} for every interior point $t$ of $I$, and \begin{align*} \Psi(t_0)=I_n. \end{align*} For each $j \in \{1,\dots,n\}$, define \begin{align*} v_j: I \to \mathbb{R}^n \end{align*} to be the $j$-th column of $\Psi$. Then $v_j$ satisfies \begin{align*} \dot{v}_j(t)=A(t)v_j(t), \qquad v_j(t_0)=e_j. \end{align*} By uniqueness in the first step, $v_j=u_j$ for every $j$. Therefore all columns of $\Psi$ and $\Phi$ agree, so $\Psi=\Phi$. [/step]

[step:Identify the solution with arbitrary initial value as $\Phi(t)x_0$] Fix $x_0 \in \mathbb{R}^n$. Define \begin{align*} x: I \to \mathbb{R}^n, \qquad t \mapsto \Phi(t)x_0. \end{align*} Since $\Phi$ is continuous on $I$ and differentiable on the interior of $I$, the map $x$ is continuous on $I$ and differentiable on the interior of $I$. For every interior point $t$ of $I$, \begin{align*} \dot{x}(t)=\Phi'(t)x_0=A(t)\Phi(t)x_0=A(t)x(t). \end{align*} Also, \begin{align*} x(t_0)=\Phi(t_0)x_0=I_n x_0=x_0. \end{align*} Thus $x(t)=\Phi(t)x_0$ is a solution of the initial value problem. If \begin{align*} y: I \to \mathbb{R}^n \end{align*} is any other solution with $y(t_0)=x_0$, then the uniqueness result in the first step, applied with $b=x_0$, gives $y=x$. Hence the unique solution with initial value $x_0$ is \begin{align*} x(t)=\Phi(t)x_0 \end{align*} for every $t \in I$. [/step]