Bounds for metric in normal coordinate Let $M$ be a Riemannian $n$-manifold and $x \in M$. For the metric tensor $g_{ij}$ of geodesic normal coordinates at $x$, there is a formula $g_{ij}(y) = \delta_{ij} + \frac13 R_{kijl} y^k y^l + O(\|y\|^3)$. Assume that the injectivity radius of $M$ is positive.
Question. Are there constants $c_1, c_2>0$ (that only depends on $M$) such that the following holds, regardless the choice of $x$ at which the geodesic normal coordinates are taken?
$$ \|y\| < c_1 \implies  | g_{ij}(y) - \delta_{ij} | < c_2 \cdot \|y\| $$
Due to higher order terms, $C$ would have to be larger than $\frac13 R_{kijl}$, but what exactly should it be? I'd like to express it in terms of simple intrinsic quantities of $M$.
The full Taylor expansion of $g_{ij}$ needs computer algebra if written in terms of derivatives of the Riemann tensor, see Equation 11.1 in here.
Related: Riemann's formula for the metric in a normal neighborhood
 A: As mentioned by Deane Yang in the comments and his (deleted) answer, one can estimate the components of the metric in normal coordinates using a transport ODE (I know it from Dolgov-Khriplovich (1983) Eq.(34), and also in more implicit form from Florides-Synge (1971), but perhaps there is a more canonical reference). The way I see the question, it does not ask to estimate the derivatives of the metric, only its pointwise deviation from the Euclidean metric $\delta_{ij}$. So I see this strategy as completely viable.
The radial ODE satisfied by the metric $g_{ij}(y)$ in normal coordinates is
$$\begin{aligned}
  r \partial_r (r\partial_r +1) (g(ry) - \delta)_{ij}
  &= \frac{1}{2} g^{kl}(ry) [r\partial_r(g(ry) - \delta)_{ik}] [r\partial_r (g(ry)-\delta)_{jl}]
    - 2r^2 y^k y^l R_{ikjl}(ry) \\
  &= \frac{1}{2} g^{kl}(ry) [r\partial_r g_{ik}(ry)] [r\partial_r g_{jl}(ry)]
    - 2r^2 y^k y^l R_{ikjl}(ry)
    =: r^2 F_{ij}(r, g(ry)-\delta, \partial_r g(ry)) ,
\end{aligned}$$
where $y=(y^i)$ is any fixed coordinate vector. The linear differential operator has the following forward Green function/fundamental solution:
$$
  G(s,r) = \Theta(s-r) \frac{1}{r} \left(1-\frac{r}{s}\right) .
$$
To solve the equation with initial condition $g_{ij}(0) = \delta_{ij}$, we can adapt the usual argument using Picard iteration and the Banach fixed point theorem to the integral equation
$$\begin{aligned}
  (g(sy)-\delta)
    &= \int_0^s dr\, r\left(1-\frac{r}{s}\right) \,  F(r, g(ry)-\delta, \partial_r g(ry))
    =: \Gamma[g(ry)-\delta] , \\
  \partial_s g(sy)
    &= \int_0^s dr \, \frac{r^2}{s^2} F(r,g(ry)-\delta,\partial_r g(ry)) ,
\end{aligned}$$
where the second equation is a consequence of the first, which we'll need shortly. Some uniform estimates on $|F|$, coupled with the integrals $\int_0^s dr\, r(1-r/s) = s^2/6$ and $\int_0^s dr\, r^2/s^2 = s/3$, will allow us to find lower bounds on the interval of existence of the solution to the transport ODE.
The original question didn't specify what norm $|g(y) - \delta|$ refers to. I'll just suppose that it is the Frobenius norm with respect to $\delta_{ij}$, $|h| = (h_{ik} h_{jl} \delta^{ij}\delta^{kl})^{1/2}$, which is quite natural for dealing with tensors in normal coordinates. To get anywhere with estimating $|F|$, it helps to assume some a priori uniform bounds on $|g(ry) - \delta| \le D < 1$, $|\partial_r g(ry)| \le C$. The constants $C$, $D$ are for now arbitrary, but they can be tuned to specific values later. The first useful uniform estimate is
$$
 |g^{kl}(ry) [r\partial_r g_{ik}(ry)] [r\partial_r g_{jl}(ry)]| \le \frac{C^2}{1-D} .
$$
Next, the curvature term satisfies
$$
  |y^k y^l R_{.k.l}(ry)| \le |R(ry)| \|y\|^2 ,
$$
where $\|y\|^2 = |yy| = y^k y^l \delta_{kl}$ and
$$
  |R(ry)| = (R_{ijkl}(ry) R_{i'k'j'l'}(ry) \delta^{ii'} \delta^{kk'} \delta^{jj'} \delta^{ll'})^{1/2}
$$ is the Frobenius norm of the endomorphism $h^{kl} \mapsto h^{kl} R_{ikjl}(ry)$. It would be nice to express the bound on the Riemann tensor in terms of some geometric quantity, but $|R(ry)|$ only has meaning inside the normal coordinate chart. However, replacing each $\delta^{ij}$ by $g^{ij}(ry)$ we get an invariant curvature scalar and we can essentially bound one with the other. Writing $\delta^{ij} = g^{ij}(ry) - (g^{ij}(ry) - \delta^{ij})$, noting that $|g^{..}(ry)-\delta^{..}| \le \frac{D}{1-D}$, and repeatedly using the Cauchy-Schwarz inequality, we get
$$\begin{aligned}
  |R(ry)| &\le \frac{\|R\|(ry) \|y\|^2}{(1-D)^2} \le \frac{\mathcal{R} \|y\|^2}{(1-D)^2} , \\
  \|R\| &= \max\{R^{(1)}, R^{(2)}, R^{(3)}, R^{(4)}\} , \\
  R^{(1)} &= (R_{ijkl} R^{ijkl})^{1/2} , \\
  R^{(2)} &= (R_{ijkl} R^{i'jkl} R_{i'j'k'l'} R^{ij'k'l'})^{1/4} , \\
  R^{(3)} &= (R_{ijkl} R^{kli'j'} R_{i'j'k'l'} R^{k'l'ij})^{1/4} , \\
  R^{(4)} &= (R_{ikjl} R^{ki'lj'} R_{i'k'j'l'} R^{k'il'j})^{1/4} , \\
  \mathcal{R} &= \sup_{x\in M} \|R\|(x) .
\end{aligned}$$
The global geometric invariant $\mathcal{R}$ of $(M,g)$ is what we will eventually use to determine the constants sought in the original question.
The main upper bound on the size of the interval $s\in [0,t]$ is that the iteration must map $(g(ry)-\delta) \mapsto \Gamma(g(ry)-\delta)$ to the space satisfying the same a priori uniform bounds. Using the above estimates on $|F|$ these inequalities are expressed as
$$\begin{aligned}
  \frac{t^2}{6} \left(\frac{1}{2} \frac{C^2}{1-D} + 2\frac{\mathcal{R \|y\|^2}}{(1-D)^2}\right) &\le D , \\
  \frac{t}{3} \left(\frac{1}{2} \frac{C^2}{1-D} + 2\frac{\mathcal{R} \|y\|^2}{(1-D)^2}\right) &\le C ,
\end{aligned}$$
which has a consistent solution $D=1/4$, $C=(2\mathcal{R}\|y\|^2/3)^{1/2}$, $t=(3\mathcal{R}\|y\|^2/8)^{-1/2}$ (no claim about the optimality of $t$). So, the integral equation $(g(ry)-\delta) = \Gamma[g(ry)-\delta]$ and the above parameters for the a priori estimates gives use the uniform bound
$$
  |g(ry)-\delta| \le \frac{r^2}{6} \left(\frac{1}{2} \frac{C^2}{1-D} + 2\frac{\mathcal{R} \|y\|^2}{(1-D)^2}\right) \le \frac{2}{3} \mathcal{R} \|ry\|^2
$$
for $\|ry\| \le (3\mathcal{R}/8)^{-1/2}$. This section on Wikipedia explains how to go from here prove the contraction property of a sufficiently high power of $\Gamma$ to prove existence.
Since now we can just replace the normal coordinate vector $ry^i$ by $y^i$ we get what was hopefully the desired estimate valid for any normal coordinate chart on $M$:
$$
  \|y\| \le \min\left\{\sqrt{\frac{8}{3\mathcal{R}}}, \rho_I\right\}
    \implies
  |g(y) - \delta| \le \frac{2}{3} \mathcal{R} \|y\|^2 ,
$$
where $\rho_I$ is the injectivity radius of $(M,g)$ (positive by assumption), which has to be there because some normal coordinate charts cannot exceed the radius $\rho_I$ anyway.
As expected, we get a larger coefficient $2/3 > 1/3$ than in the Taylor formula. It might be possible to optimize this coefficient and the bound on $\|y\|$ a bit more, if needed. Also changing the Frobenius norm to a different one could also affect the constant and the definition of $\mathcal{R}$ in terms of global geometric invariants of $(M,g)$.
A: APOLOGIES: Major revision, including details, below. My original answer was incorrect.
Thanks to @IgorKhavkine for pointing out that the question asks for a point wise bound on only the metric itself and not its first derivative, as well as his impressive analysis of the ODE satisfied by the metric written in geodesic normal coordinates. I'd like to describe a different approach using Jacobi fields.
Start with an orthonormal basis $(\partial_1, \dots, \partial_n)$ of $T_xM$. For each $v = v^i\partial_i \in T_xM$,
$$
\exp_x(v) = \gamma(1),
$$
where $\gamma$ is the constant speed geodesic satisfying $\gamma(0) = x$ and $\gamma'(0) = v$.
If we denote $\Phi(v^1, \dots, v^n) = \exp_x(v)$, then the metric in geodesic normal coordinates is given by
$$ g_{ij} = g(\partial_i\Phi,\partial_j\Phi). $$
On the other hand, if we let $\gamma_t$ be the constant speed geodesic with $\gamma_t(0) = x$ and $\gamma_t'(0) = v + t\partial_i$, then
$$
\partial_i\Phi(v) = \left.\frac{d}{dt}\right|_{t=0}\Phi(v+te_i) = \left.\frac{d}{dt}\right|_{t=0})\gamma_t(1) = J_i(1), $$
where
$$
J_i(s) = \left.\frac{d}{dt}\right|_{t=0}\gamma_t(s)$$
is the unique Jacobi field that satisfies
$$ J_i(0) = 0\text{ and } \nabla_{v}J_i(0) = \partial_i.
$$
Define an orthonormal frame $(e_1, e_n)$ by parallel translating $(\partial_1, \dots, \partial_n)$ along each geodesic passing through $x$ and, abusing notation, let $J$ denote the matrix whose components $J_i^j$ are given by
$$ J_i = J_i^je_j. $$
The Jacobi equation,
$$\nabla_T\nabla TJ_i = R(T,J_i)T,\ 1 \le i \le n,$$
where $T = \gamma'(t)$, is equivalent to
$$
J'' + JK = 0,
$$
where $K$ is the matrix such that (up to sign) $$K^i_j = e_i\cdot R(e_j,T)T.$$
We therefore want to get pointwise bounds on $J(1) \in \mathbb{R}^n$ satisfying the ODE above with $J(0) = 0$ and $J'(0) = I$. This is now straightforward. Integrating the ODE twice, we get
$$ J(t) = tI - \int_{s=0}^{s=t} (t-s)JK\,ds. $$
In particular,
$$ \partial \Phi(v) = J(1) = I - \int_{s=0}^{s=1} (1-s)JK\,ds. $$
The desired bound now follows easily.
