Assume we have a real-valued function $f=f_\epsilon(x)$ depending on a (typically small) parameter $\epsilon$, whose inverse (with respect to $x$) is well-defined but hard to compute. What can be said about $\frac{\partial}{\partial \epsilon}f^{-1}$ (without computing $f^{-1}$ explicitly)?
For the sake of concreteness, consider $f_\epsilon(x):=x\sqrt{1+\epsilon^2(1+x+x^2)}.$ Then f'_\epsilon(x)>0 for every $x\in\mathbb{R}$, and so $g_\epsilon:=f_\epsilon^{-1}$ is well defined as a function of $x$ on the whole real line $\mathbb{R}$. An exact expression for $f_\epsilon^{-1}$ would involve solving a quartic equation in $x$, which we don't want to do unless it turns out to be really necessary. My question is: how can we get some control on $\frac{\partial}{\partial \epsilon}g_\epsilon$? More precisely, it is straightforward to verify that $\frac{\partial}{\partial \epsilon}f_\epsilon(x)=\frac{\epsilon x (1+x+x^2)}{\sqrt{1+\epsilon^2(1+x+x^2)}},$ and so $\Big|\frac{\partial}{\partial \epsilon}f_\epsilon(x)\Big|\leq C_1(1+x+x^2)$ as long as $x$ lies in the bounded region $|x|\leq C_1\epsilon^{-1}$. These are the kind of (uniform in $\epsilon$) bounds I am interested in. Can one find similar (uniform in $\epsilon$) bounds for $\Big|\frac{\partial}{\partial \epsilon}g_\epsilon(y)\Big|$ which hold in a region $|y|\leq C_2\epsilon^{-1}$, say?
UPDATE: I was making the mistake of thinking that one could use the chain rule to evaluate
\begin{equation} \textrm{(1)}\quad\quad\quad \frac{\partial}{\partial\epsilon}(F_\epsilon\circ G_\epsilon)(x) \end{equation} (here $F_\epsilon$ and $G_\epsilon$ are just some functions of $x$ depending on $\epsilon$, not necessarily inverses of each other) as a product of two terms. That turns out not to be the case, but one should still be able to express (1) as a sum of two or more terms. Any ideas?
Thank you.