In Giaquinta's and Hildebrandt's 1996, "Calculus of Variations 1", pages 147-148, they develop the definition of inner variations.
They first fix $\lambda\in C_c^{\infty}(\overline\Omega;\mathbb{R}^n)$ and $0<\varepsilon_0<(2K)^{-1}$, where $\Omega\subset\mathbb{R}^n$ is an open bounded set and $K$ is the Lipschitz constant of $\lambda$.
Then they show that the family of maps: \begin{equation} \mathcal{F}=\{\xi_{\varepsilon}: \overline\Omega\rightarrow\overline\Omega\ |\ \varepsilon\in (-\varepsilon_0, \varepsilon_0)\} \end{equation} where \begin{equation} \xi_{\varepsilon}(y)=y+\varepsilon \lambda(y)\quad y\in \overline\Omega, \ \varepsilon\in(-\varepsilon_0, \varepsilon_0) \end{equation} is actually a family of diffeomorphisms of $\overline\Omega$ that further satisfy the conditions:
- $\xi_0(y)=y$,
- $\frac{\partial}{\partial\varepsilon}\xi_{0}(y)=\lambda(y)$ and
- $\xi_{\varepsilon}(y)=y\quad \forall\ x\in\partial\Omega, |\varepsilon|<\varepsilon_0$.
Now, I follow their arguments for demonstrating the above. However, they further say that the inverse mapping of $\xi_{\varepsilon}$ is given by: \begin{equation*} \eta_{\varepsilon}(x)=x-\varepsilon\lambda(x)+o(\varepsilon) \end{equation*}but I don't understand this because if $y\in\text{spt}(\lambda)$ then one would expect that $\eta_{\varepsilon}\circ\xi_{\varepsilon}(y)=y$. However, I get as far as:
\begin{equation*} \eta_{\varepsilon}\circ\xi_{\varepsilon}(y)=y+\varepsilon\lambda(y)-\varepsilon\lambda[y+\varepsilon\lambda(y)]+o(\varepsilon). \end{equation*}
Since $\lambda$ is continuous, we have $\lambda(y+\varepsilon \lambda(y)) = \lambda(y)+o(1)$. Here $o(1)$ (as $\varepsilon\to 0$) is uniform in $y$, since $\lambda$ is a compactly supported perturbation. Hence, \begin{equation*} y+\varepsilon\lambda(y)-\varepsilon\lambda[y+\varepsilon\lambda(y)]+o(\varepsilon) =y+\varepsilon\lambda(y)-\varepsilon\lambda(y) + o( \varepsilon) +o(\varepsilon) =y+o(\varepsilon) \end{equation*} as expected.
Another way to approach inner variation is to take $\xi$ to be the flow map associated to vector field $\lambda$. Then $\xi_\varepsilon\circ \xi_{-\varepsilon}$ is obviously the identity map, and $\xi_\varepsilon(y)=y+\varepsilon \lambda(y)+o(\varepsilon)$ follows from the continuity of $\lambda$. This way, positive and negative values of $\varepsilon$ are treated simultaneously.
You may want to consult Two-dimensional geometric variational problems by Jost, which also has a derivation of inner variational equations (though in two dimensions only).