[proofplan]
The proof reduces the first-order PDE to an ordinary differential equation along the characteristic curve. After defining $v(t)=u(X(t))$, the chain rule and the characteristic equation give $v'(t)+c(X(t))v(t)=f(X(t))$. We then multiply by the standard integrating factor based at $t_0$, integrate the resulting exact derivative, and rearrange the expression to obtain the asserted representation.
[/proofplan]
[step:Convert the PDE into a scalar equation along the characteristic]
Define the composed function
\begin{align*}
v:I \to \mathbb{R}
\end{align*}
by $v(t)=u(X(t))$ for $t \in I$. Since $u \in C^1(U)$ and $X \in C^1(I;U)$, the function $v$ belongs to $C^1(I)$. For each $t \in I$, the chain rule gives
\begin{align*}
v'(t)=\nabla u(X(t))\cdot X'(t).
\end{align*}
Because $X$ is a characteristic curve for $a$, we have $X'(t)=a(X(t))$. Therefore
\begin{align*}
v'(t)=\nabla u(X(t))\cdot a(X(t)).
\end{align*}
Evaluating the PDE at the point $X(t)\in U$ gives
\begin{align*}
a(X(t))\cdot \nabla u(X(t))+c(X(t))u(X(t))=f(X(t)).
\end{align*}
Since the Euclidean dot product is symmetric and $v(t)=u(X(t))$, this becomes
\begin{align*}
v'(t)+c(X(t))v(t)=f(X(t)).
\end{align*}
Thus $v$ satisfies the inhomogeneous scalar linear equation along the curve.
[/step]
[step:Multiply by the integrating factor based at $t_0$]
Define the coefficient along the characteristic
\begin{align*}
\beta:I \to \mathbb{R}
\end{align*}
by $\beta(r)=c(X(r))$ for $r \in I$. Since $c \in C^0(U)$ and $X \in C^1(I;U)$, the function $\beta$ is continuous on $I$.
Fix $t_0 \in I$. Define the integrating factor
\begin{align*}
M:I \to (0,\infty)
\end{align*}
by
\begin{align*}
M(s)=\exp\left(\int_{t_0}^{s}\beta(r)\,d\mathcal{L}^1(r)\right).
\end{align*}
The function $M$ is $C^1$ and satisfies
\begin{align*}
M'(s)=\beta(s)M(s)
\end{align*}
for every $s \in I$. Multiplying the scalar equation
\begin{align*}
v'(s)+\beta(s)v(s)=f(X(s))
\end{align*}
by $M(s)$ gives
\begin{align*}
M(s)v'(s)+M(s)\beta(s)v(s)=M(s)f(X(s)).
\end{align*}
Using $M'(s)=\beta(s)M(s)$, the left-hand side is the product derivative:
\begin{align*}
\frac{d}{ds}\bigl(M(s)v(s)\bigr)=M(s)f(X(s)).
\end{align*}
[guided]
The goal is to remove the lower-order term $\beta(s)v(s)$ from the scalar equation. Define
\begin{align*}
\beta:I \to \mathbb{R}
\end{align*}
by $\beta(r)=c(X(r))$. This is continuous because $c$ is continuous on $U$ and $X$ is continuous as a $C^1$ map into $U$.
We now choose a positive function $M$ whose derivative produces exactly the coefficient multiplying $v$. Fix $t_0 \in I$ and define
\begin{align*}
M:I \to (0,\infty)
\end{align*}
by
\begin{align*}
M(s)=\exp\left(\int_{t_0}^{s}\beta(r)\,d\mathcal{L}^1(r)\right).
\end{align*}
Since $\beta$ is continuous, the function $s \mapsto \int_{t_0}^{s}\beta(r)\,d\mathcal{L}^1(r)$ is $C^1$ with derivative $\beta(s)$. Differentiating the exponential composition gives
\begin{align*}
M'(s)=\beta(s)M(s).
\end{align*}
Along the characteristic, the previous step proved
\begin{align*}
v'(s)+\beta(s)v(s)=f(X(s)).
\end{align*}
Multiplying this identity by $M(s)$ gives
\begin{align*}
M(s)v'(s)+M(s)\beta(s)v(s)=M(s)f(X(s)).
\end{align*}
The reason for the choice of $M$ is now visible: because $M'(s)=\beta(s)M(s)$, the left-hand side is exactly the product derivative of $M(s)v(s)$. Hence
\begin{align*}
\frac{d}{ds}\bigl(M(s)v(s)\bigr)=M(s)f(X(s)).
\end{align*}
This turns the inhomogeneous linear equation into an equation that can be integrated directly.
[/guided]
[/step]
[step:Integrate the exact derivative and solve for $v(t)$]
Let $t \in I$ be arbitrary. Integrating the identity
\begin{align*}
\frac{d}{ds}\bigl(M(s)v(s)\bigr)=M(s)f(X(s))
\end{align*}
from $t_0$ to $t$ gives
\begin{align*}
M(t)v(t)-M(t_0)v(t_0)=\int_{t_0}^{t}M(s)f(X(s))\,d\mathcal{L}^1(s).
\end{align*}
By definition of $M$, we have $M(t_0)=1$. Therefore
\begin{align*}
M(t)v(t)=v(t_0)+\int_{t_0}^{t}M(s)f(X(s))\,d\mathcal{L}^1(s).
\end{align*}
Since $M(t)>0$, divide by $M(t)$:
\begin{align*}
v(t)=v(t_0)M(t)^{-1}+\int_{t_0}^{t}M(t)^{-1}M(s)f(X(s))\,d\mathcal{L}^1(s).
\end{align*}
Now
\begin{align*}
M(t)^{-1}=\exp\left(-\int_{t_0}^{t}\beta(r)\,d\mathcal{L}^1(r)\right),
\end{align*}
and, for each $s \in I$,
\begin{align*}
M(t)^{-1}M(s)=\exp\left(\int_{t_0}^{s}\beta(r)\,d\mathcal{L}^1(r)-\int_{t_0}^{t}\beta(r)\,d\mathcal{L}^1(r)\right).
\end{align*}
By additivity of oriented integrals,
\begin{align*}
\int_{t_0}^{s}\beta(r)\,d\mathcal{L}^1(r)-\int_{t_0}^{t}\beta(r)\,d\mathcal{L}^1(r)=-\int_s^t\beta(r)\,d\mathcal{L}^1(r).
\end{align*}
Thus
\begin{align*}
M(t)^{-1}M(s)=\exp\left(-\int_s^t\beta(r)\,d\mathcal{L}^1(r)\right).
\end{align*}
Substituting $\beta(r)=c(X(r))$, $v(t)=u(X(t))$, and $v(t_0)=u(X(t_0))$ yields
\begin{align*}
u(X(t))=u(X(t_0))\exp\left(-\int_{t_0}^{t}c(X(r))\,d\mathcal{L}^1(r)\right)+\int_{t_0}^{t}\exp\left(-\int_s^t c(X(r))\,d\mathcal{L}^1(r)\right)f(X(s))\,d\mathcal{L}^1(s).
\end{align*}
This is the desired inhomogeneous [representation formula](/theorems/39) along the characteristic.
[/step]