As Thomas Rot already suggested: this follows directly from smooth dependence of ODEs on initial conditions.

Let $p \in M$ and denote by $\Phi^t\colon TM \to TM$ the geodesic flow. Then the exponential map at $p \in M$ is defined as the time one geodesic flow, restricted to $T_p M$ and projected onto $M$, i.e.
$$
  \exp_p = \pi \circ \Phi^1|_{T_p M} \colon T_p M \to M
$$
where $\pi$ is the tangent bundle projection.

In local coordinates around $p$ this amounts to an ODE on $\mathbb{R}^{2n}$ involving the Christoffel symbols, and these are smooth since $(M,g)$ was assumed smooth.

**Edit added:** In local coordinates the geodesic flow is given by
$$
  \dot{x}^i = v^i, \qquad \dot{v}^i = -\Gamma^i_{jk}(x) v^j v^k
$$
with $(x,v) \in \mathbb{R}^{2n}$ induced local coordinates on $T M$.
This is well-defined, also for $v = 0$.