The most general form of a linear IVP that was considered in my course is $$\dot x(t) = A(t)x(t)+b(t),\quad t\in J, \quad x(t_0)=x_0,$$ for $J$ an interval, $t_0\in J$, $A\in C(J,\mathbb{R}^{m\times m})$, and $b\in C(J,\mathbb{R}^m)$. The unique solution is derived using fundamental matrices and given as $$x(t) = X(t)\left(X^{-1}(t_0)x_0+\int_{t_0}^t X^{-1}(\tau)b(\tau)\operatorname{d}\tau\right)$$ with $X$ being a fundamental matrix.
However, considering the same ODE in general Banach spaces, that is, $A\in C(J, \mathcal{L}(E))$, $b\in C(J, E)$, $E$ being a banach-space, what's the change? Is $$(\star)\quad x(t) = e^{\int_{t_0}^t A(\tau)\operatorname{d}\tau}x_0 +\int_{t_0}^t e^{\int_{\tau}^t A(s)\operatorname{d} s}b(\tau)\operatorname{d}\tau$$ the unique solution to the IVP? I understand the possibility of the space of solutions being of infinite dimension prohibits the use of fundamental matrices. Are fundamental matrices used to give concrete answers to concrete IVP's only?
As I have worked through linear ODE's with constant $A\in\mathcal{L}(E)$ in banach-spaces, what's the difference coming from those ODE's? What becomes false after replacing "$(t-t_0)A$" by "$\int_{t_0}^tA(\tau)\operatorname{d}\tau$"?
Yes, it's still the unique solution. The same arguments that worked in the finite dimensional case still work since the contraction mapping principle is true in arbitrary complete metric spaces. (It doesn't use Heine-Borel.) We do use compactness of closed intervals in $\mathbb{R}$, though, to get the necessary Lipschitz bounds on the vector field $A$. I give a complete proof below.
Assume without loss of generality that $J \ni 0$. Suppose $K$ is a compact subinterval of $J$ containing $0$. Then the mappings $A \restriction_{K} : K \to \mathcal{L}(E)$ and $b \restriction_{K} : K \to E$ are bounded, since their ranges are compact. Consider the solution mapping $$\Phi \quad : \quad \alpha \mapsto \alpha(0) + \int_{0}^{t} \left(A(s)\alpha(s) + b(s)\right) \, ds.$$
We need to define a range for $\Phi$. We can do this as in the finite dimensional case. Suppose the initial condition we're interested in is $x_{0} \in E$ and set $R = 2 \|x_{0}\|$. Set $\ell_{1} = \frac{1}{2 \|A\|_{K}}$, where $\|A\|_{K} = \max \{\|A(t)\| \, \mid \, t \in K\}$, and set $\ell_{2} = \frac{R}{2} \cdot \left(2 R\|A\|_{K} + \|b\|_{K}\right)^{-1}$ with $\|b\|_{K} = \max\{\|b(t)\| \, \mid \, t \in K\}$. Let $\ell = \min\{\ell_{1},\ell_{2},L(K)\}$, where $L > 0$ satisfies $[-L,L] \subseteq K$. Let $\mathscr{C} = \{\alpha \in C([-\ell,\ell],E) \, \mid \, \alpha(0) = x_{0}, \, \, \|\alpha\|_{C([-\ell,\ell],E)} \leq R\}$, where $\|\alpha\|_{C([-\ell,\ell],E)} = \max\{\|\alpha(t)\| \, \mid \, t \in [-\ell,\ell]\}$.
We claim $\Phi : \mathscr{C} \to \mathscr{C}$ is a well-defined contraction mapping. Since $\mathscr{C}$ is a closed subset of a Banach space, it's a complete metric space. Therefore, the contraction mapping principle implies there is a unique fixed point $\tilde{\alpha} \in \mathscr{C}$, which satisfies $$\forall t \in [-\ell,\ell] \quad \tilde{\alpha}(t) = \tilde{\alpha}(0) + \int_{0}^{t} (A(s) \tilde{\alpha}(s) + b(s)) \, ds.$$
Since the integrand is continuous, $\tilde{\alpha} \in C^{1}([-\ell,\ell],E)$ and $\frac{d \tilde{\alpha}}{d t} = A(t) \tilde{\alpha}(t) + b(t)$. By uniqueness, the Duhamel formula you provided is the solution (at least on $[-\ell,\ell]$). It remains to prove $\mathscr{C}$ is well-defined and a contraction mapping.
First, we show it's well-defined. If $t \in [-\ell,\ell]$ and $\alpha \in \mathscr{C}$, then \begin{align*} \left\|\alpha(0) + \int_{0}^{t} \left(A(s) \alpha(s) + b(s)\right) \, ds \right\| &\leq \|\alpha(0)\| + \|A\|_{K} \int_{0}^{t} \|\alpha(s)\| \, ds + \int_{0}^{t} \|b(s)\| \, ds \\ &\leq \frac{R}{2} + \left(\|A\|_{K} R + \|b\|_{K}\right)\ell \\ &\leq \frac{R}{2} + \left(\|A\|_{K} R + \|b\|_{K}\right) \ell_{2} \\ &\leq R. \end{align*} This shows $\Phi$ is a well-defined operator on $\mathscr{C}$.
Next, to see that $\Phi$ is a contraction mapping, observe that, if $\alpha,\beta \in \mathscr{C}$ and $t \in [-\ell,\ell]$, then \begin{align*} \|\alpha(t) - \beta(t)\| &= \left \|\int_{0}^{t} A(s)(\alpha(s) - \beta(s)) \, ds \right \| \\ &\leq \|A\|_{K} \int_{0}^{t} \|\alpha(s) - \beta(s)\| \, ds \\ &\leq \|A\|_{K} \ell \|\alpha - \beta\|_{C([-\ell,\ell],E)} \\ &\leq \frac{1}{2} \|\alpha - \beta\|_{C([-\ell,\ell],E)}. \end{align*} Thus, $\|\Phi(\alpha) - \Phi(\beta)\|_{C([-\ell,\ell],E)} \leq \frac{1}{2} \|\alpha - \beta\|_{C([-\ell,\ell],E)}$, proving the claim.
The above proves we can find a unique solution on $[-\ell,\ell]$. However, the choice of $\ell$ works uniformly on $K$ so we can repeat the same argument to extend the solution to $[\ell,2 \ell]$, $[2 \ell, 3 \ell]$, and so on, as long as $[n \ell, (n + 1) \ell] \subseteq K$. Thus, we will be able to solve the ODE uniquely on $K$. Since $J$ equals the union of its compact subintervals, this gives uniqueness on the original interval $J$.
A question I don't know the answer to but there is literature on is the case of the ODE $\dot{x} = Ax$, where $A$ is an unbounded linear operator. It might be worth looking into this if you have some familiarity with functional analysis since partial differential equations can be studied from this viewpoint.