[proofplan]
We first rewrite both solutions in integral form on the interval $[a,b]$ and subtract the two identities. The assumed Lipschitz bound on $[a,b]\times C$ then converts the difference equation into a scalar integral inequality for $u(t)=|x(t)-y(t)|$. Finally, we prove the needed Gronwall estimate directly by applying the product rule to an exponentially weighted auxiliary function.
[/proofplan]
[step:Subtract the integral identities for the two solutions]
Throughout the proof, $\mathcal{L}^1$ denotes one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}$.
Define the difference map
\begin{align*}
d:[a,b]\to \mathbb{R}^n,\qquad d(t)=x(t)-y(t).
\end{align*}
Since $F$ is continuous, $[a,b]\subset I$, and $x$ and $y$ are solutions of $\dot z=F(t,z)$ on $[a,b]$, the integral formulation of solutions from [citetheorem:8311] gives, for every $t\in [a,b]$,
\begin{align*}
x(t)=x(a)+\int_a^t F(s,x(s))\,d\mathcal{L}^1(s)
\end{align*}
and
\begin{align*}
y(t)=y(a)+\int_a^t F(s,y(s))\,d\mathcal{L}^1(s).
\end{align*}
Subtracting these two vector identities gives
\begin{align*}
d(t)=d(a)+\int_a^t \bigl(F(s,x(s))-F(s,y(s))\bigr)\,d\mathcal{L}^1(s).
\end{align*}
[guided]
We want to compare $x$ and $y$ at time $t$ using only their initial separation at time $a$. The right way to expose that dependence is to write each solution as its starting value plus accumulated velocity.
Throughout the proof, $\mathcal{L}^1$ denotes one-dimensional Lebesgue measure on $\mathbb{R}$. Define
\begin{align*}
d:[a,b]\to \mathbb{R}^n,\qquad d(t)=x(t)-y(t).
\end{align*}
The maps $x$ and $y$ are solutions of $\dot z=F(t,z)$ on $[a,b]$. Therefore the integral formulation of an initial value problem, [citetheorem:8311], applies to each of them with initial time $a$. For every $t\in [a,b]$, it gives
\begin{align*}
x(t)=x(a)+\int_a^t F(s,x(s))\,d\mathcal{L}^1(s)
\end{align*}
and
\begin{align*}
y(t)=y(a)+\int_a^t F(s,y(s))\,d\mathcal{L}^1(s).
\end{align*}
The hypotheses needed here are exactly the continuity of $F$, the fact that $x$ and $y$ solve the differential equation on the interval, and the inclusion $[a,b]\subset I$.
Subtracting the second identity from the first, using linearity of the vector-valued integral, yields
\begin{align*}
d(t)=d(a)+\int_a^t \bigl(F(s,x(s))-F(s,y(s))\bigr)\,d\mathcal{L}^1(s).
\end{align*}
This identity is the starting point for the estimate: it says that the later separation equals the initial separation plus the accumulated difference of the vector field along the two trajectories.
[/guided]
[/step]
[step:Convert the difference identity into a scalar integral inequality]
Define
\begin{align*}
u:[a,b]\to [0,\infty),\qquad u(t)=|d(t)|=|x(t)-y(t)|.
\end{align*}
The map $u$ is continuous because $x$ and $y$ are continuous. Taking Euclidean norms in the identity for $d(t)$, applying the triangle inequality for vector-valued integrals, and using $x(s),y(s)\in C$ for every $s\in[a,b]$, we obtain
\begin{align*}
u(t)\le u(a)+\int_a^t |F(s,x(s))-F(s,y(s))|\,d\mathcal{L}^1(s).
\end{align*}
Since $L$ is a Lipschitz constant for $F$ in the state variable on $[a,b]\times C$,
\begin{align*}
|F(s,x(s))-F(s,y(s))|\le L|x(s)-y(s)|=Lu(s)
\end{align*}
for every $s\in[a,b]$. Hence, for every $t\in[a,b]$,
\begin{align*}
u(t)\le u(a)+L\int_a^t u(s)\,d\mathcal{L}^1(s).
\end{align*}
[/step]
[step:Prove the scalar Gronwall estimate with an exponential weight]
Let
\begin{align*}
\alpha=u(a)=|x(a)-y(a)|.
\end{align*}
Fix $\varepsilon>0$ and define
\begin{align*}
h_\varepsilon:[a,b]\to (0,\infty),\qquad h_\varepsilon(t)=\alpha+\varepsilon+L\int_a^t u(s)\,d\mathcal{L}^1(s).
\end{align*}
The function $h_\varepsilon$ is continuously differentiable on $(a,b)$, and for every $t\in(a,b)$,
\begin{align*}
h_\varepsilon'(t)=Lu(t).
\end{align*}
The integral inequality gives $u(t)\le h_\varepsilon(t)$ for every $t\in[a,b]$, so
\begin{align*}
h_\varepsilon'(t)\le Lh_\varepsilon(t)
\end{align*}
for every $t\in(a,b)$. Define
\begin{align*}
r_\varepsilon:[a,b]\to (0,\infty),\qquad r_\varepsilon(t)=e^{-L(t-a)}h_\varepsilon(t).
\end{align*}
For $t\in(a,b)$, the product rule gives
\begin{align*}
r_\varepsilon'(t)=e^{-L(t-a)}\bigl(h_\varepsilon'(t)-Lh_\varepsilon(t)\bigr)\le 0.
\end{align*}
Thus $r_\varepsilon$ is nonincreasing on $[a,b]$, and therefore
\begin{align*}
e^{-L(t-a)}h_\varepsilon(t)=r_\varepsilon(t)\le r_\varepsilon(a)=\alpha+\varepsilon
\end{align*}
for every $t\in[a,b]$. Multiplying by $e^{L(t-a)}$ gives
\begin{align*}
h_\varepsilon(t)\le e^{L(t-a)}(\alpha+\varepsilon).
\end{align*}
Since $u(t)\le h_\varepsilon(t)$, we have
\begin{align*}
u(t)\le e^{L(t-a)}(\alpha+\varepsilon).
\end{align*}
Letting $\varepsilon\downarrow 0$ gives
\begin{align*}
u(t)\le e^{L(t-a)}\alpha.
\end{align*}
[guided]
We now prove the precise scalar estimate needed from the inequality
\begin{align*}
u(t)\le \alpha+L\int_a^t u(s)\,d\mathcal{L}^1(s),
\end{align*}
where
\begin{align*}
\alpha=u(a)=|x(a)-y(a)|.
\end{align*}
The small parameter $\varepsilon>0$ avoids any division or positivity issue when $\alpha=0$.
Fix $\varepsilon>0$ and define the auxiliary function
\begin{align*}
h_\varepsilon:[a,b]\to (0,\infty),\qquad h_\varepsilon(t)=\alpha+\varepsilon+L\int_a^t u(s)\,d\mathcal{L}^1(s).
\end{align*}
Because $u$ is continuous, the one-variable [fundamental theorem of calculus](/theorems/632) implies that $h_\varepsilon$ is continuously differentiable on $(a,b)$ and satisfies
\begin{align*}
h_\varepsilon'(t)=Lu(t)
\end{align*}
for every $t\in(a,b)$. The scalar integral inequality says
\begin{align*}
u(t)\le \alpha+L\int_a^t u(s)\,d\mathcal{L}^1(s)\le h_\varepsilon(t).
\end{align*}
Therefore
\begin{align*}
h_\varepsilon'(t)=Lu(t)\le Lh_\varepsilon(t)
\end{align*}
for every $t\in(a,b)$.
The exponential factor is chosen so that differentiating cancels the term $Lh_\varepsilon(t)$. Define
\begin{align*}
r_\varepsilon:[a,b]\to (0,\infty),\qquad r_\varepsilon(t)=e^{-L(t-a)}h_\varepsilon(t).
\end{align*}
Using the product rule, for every $t\in(a,b)$,
\begin{align*}
r_\varepsilon'(t)=e^{-L(t-a)}h_\varepsilon'(t)-Le^{-L(t-a)}h_\varepsilon(t).
\end{align*}
Factoring the common positive term gives
\begin{align*}
r_\varepsilon'(t)=e^{-L(t-a)}\bigl(h_\varepsilon'(t)-Lh_\varepsilon(t)\bigr)\le 0.
\end{align*}
Hence $r_\varepsilon$ is nonincreasing on $[a,b]$. Consequently, for every $t\in[a,b]$,
\begin{align*}
e^{-L(t-a)}h_\varepsilon(t)=r_\varepsilon(t)\le r_\varepsilon(a)=h_\varepsilon(a)=\alpha+\varepsilon.
\end{align*}
Multiplying by the positive number $e^{L(t-a)}$ yields
\begin{align*}
h_\varepsilon(t)\le e^{L(t-a)}(\alpha+\varepsilon).
\end{align*}
Since $u(t)\le h_\varepsilon(t)$, we obtain
\begin{align*}
u(t)\le e^{L(t-a)}(\alpha+\varepsilon).
\end{align*}
This estimate holds for every $\varepsilon>0$, so taking the limit as $\varepsilon\downarrow 0$ gives
\begin{align*}
u(t)\le e^{L(t-a)}\alpha.
\end{align*}
[/guided]
[/step]
[step:Translate the scalar estimate back to the two solutions]
By the definitions of $u$ and $\alpha$, the preceding estimate says that for every $t\in[a,b]$,
\begin{align*}
|x(t)-y(t)|\le e^{L(t-a)}|x(a)-y(a)|.
\end{align*}
This is the desired stability estimate for the two solutions on $[a,b]$.
[/step]