[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$.
[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]
[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]