I have some doubts about this discretization of the geodesic flow. I believe you should have another component of the map: one that transforms the velocity vector $V$. By the way, technically $V$ should be a momentum covector (belonging to the cotangent bundle).
If you want a discrete version of the geodesic flow, I think you should start from a Lagrangian point of view. The geodesics on a Riemannian manifold with metric tensor $G(x)$ are the solutions of the Euler-Lagrange equations $\frac{d}{dt}\big(\frac{\partial L}{\partial \dot{x}}(x,\dot{x})\big) = \frac{\partial L}{\partial x}(x,\dot{x})$ with Lagrange function $L(x,\dot{x}) = \frac{1}{2}\dot{x}^T G(x) \dot{x}$. Now, a straight-forward discretization of this method is as follows: replace the derivative $\dot{x}$ by a difference $(x_1 - x)/\varepsilon$ and consider the discrete Lagrangian $L_{\varepsilon}(x,x_1) := L(x, (x_1-x)/\varepsilon)$. By analogy with the continuous case, where the action is $S[x] = \int_{t_0}^{t_1} L(x,\dot{x})dt$ , in the discrete case the action is $S_{\varepsilon}[x] = \sum_{k=0}^{N} L_{\varepsilon}(x_k,x_{k+1})$. Then, the critical discrete geodesic should simply be the solution to the algebraic equations $\nabla S[x] = 0$ which componentwise is basically the discrete Euler-Lagrange equations $\frac{\partial L_{\varepsilon}}{\partial x_k}(x_{k-1}, x_k) + \frac{\partial L_{\varepsilon}}{\partial x_k}(x_{k}, x_{k+1}) = 0$ for $k=1,...,N$. Now, we can look locally at these discrete equations and write them down using a simplified subscript notation $\frac{\partial L_{\varepsilon}}{\partial x}(x_{-1}, x) + \frac{\partial L_{\varepsilon}}{\partial x}(x, x_{1}) = 0$. Introduce the variable $p = \frac{\partial L_{\varepsilon}}{\partial x}(x, x_{1})$. Then the discrete Euler-Lagrange equations become $\frac{\partial L_{\varepsilon}}{\partial x}(x_{-1}, x) +p= 0$ which shifted by one subscript turn into $\frac{\partial L_{\varepsilon}}{\partial x_1}(x, x_1) + p_1= 0$. Thus we obtain the equations $$p = \frac{\partial L_{\varepsilon}}{\partial x}(x, x_{1})$$
$$p_1 = -\frac{\partial L_{\varepsilon}}{\partial x_1}(x, x_1).$$ If one can express $x_1$ as a function of $(x,p)$ from the first equation, then the second equation also gives us $p_1$ as a function of $(x,p)$. Thus we have obtained a map $\Phi : (x,p) \mapsto (x_1,p_1)$. Observe that this is a map $\Phi : T^*M \to T^*M$, which turns out to be symplectic, because the Lagrangian $L_{\varepsilon}$ is in fact a generating function of the symplectomorphism $\Phi$. In the current situation $L_{\varepsilon}(x,x_1) = \frac{1}{2 \varepsilon^2} (x_1 - x)^T G(x) (x_1 - x)$, so the second equation above yields $p_1 = -\frac{1}{\varepsilon^2} G(x)(x_1-x)$, which is equivalent to $x_1 = x - \varepsilon^2 G(x)^{-1}p_1$, which resembles your equation. But as you can see, it does not involve a constant tangent vector $V$, but rather a variable cotangent vector $p_1$, which can be found by taking into account the first equation as well. And in this construction you have convergency of $x_n : n =0,1,2,...$ to a smooth geodesic.
