Taylor series (or equivalent at $\epsilon\to0$) of the integral over $x$ of a function of $x$ and $\epsilon$

This uses the same basic idea as Christian Blatter's answer. Observe that $$ f(x,\epsilon) \;=\; g_3(x)\epsilon^3 \,+\, g_4(x)\epsilon^4 \,+\, g_5(x)\epsilon^5 \,+\, \cdots $$ where each $g_n(x)$ is the product of a polynomial with $e^{-x^2/2}$: $$ \begin{align*} g_3(x) \;&=\; \frac12 \bigl(-x^3 + x^2 + 5x -1\bigr)e^{-x^2/2} \\[6pt] g_4(x) \;&=\; \frac{1}{24} \bigl(-3x^5+46x^3-12x^2-93x+12)e^{-x^2/2} \\[6pt] g_5(x) \;&=\; \frac{1}{48}\bigl(-x^7 + 27 x^5 - 10 x^4 - 165 x^3 + 84 x^2 + 195 x - 54\bigr)e^{-x^2/2} \\[6pt] g_6(x) \;&=\; \left(\tfrac{- 5 x^9 + 220 x^7 - 120 x^6 - 2714 x^5 + 2280 x^4 + 10020 x^3 - 8280 x^2 - 7725 x + 3240}{1920}\right)e^{-x^2/2}\\ &\;\vdots \end{align*} $$ The function $f(x,\epsilon)$ is zero along three curves:

These curves hit the $x$-axis at the three roots of the polynomial $x^3 - x^2 - 5x + 1$. It is easy to check that $\dfrac{\partial}{\partial \epsilon}[f(x,\epsilon)/\epsilon^3]\ne 0$ at these three points, so the picture really does look like this in some neighborhood of the $x$-axis, and the intersections really are transverse.

If we write the three curves as equations of the form $$ x \;=\; \psi_1(\epsilon),\qquad x=\psi_2(\epsilon),\qquad x=\psi_3(\epsilon) $$ Note that the functions $\psi_1$, $\psi_2$, and $\psi_3$ are $C^\infty$, by the Implicit Function Theorem. Then $$ \Delta(\epsilon) \;=\; \int_{-\infty}^{\psi_1(\epsilon)} \!\!f(x,\epsilon)\,dx \,-\, \int_{\psi_1(\epsilon)}^{\psi_2(\epsilon)} \!\!f(x,\epsilon)\,dx \,+\, \int_{\psi_2(\epsilon)}^{\psi_3(\epsilon)} \!\!f(x,\epsilon)\,dx \,-\, \int_{\psi_3(\epsilon)}^{\infty} \!\!f(x,\epsilon)\,dx $$ We can use this to derive the first few terms of a Taylor expansion for $\Delta(\epsilon)$.

Specifically, we have $$ \Delta(\epsilon) \;=\; \kappa \epsilon^3 \,+\, \lambda \epsilon^4 \,+\,\mu\epsilon^5 \,+\, o(\epsilon^5) $$ for some constants $\kappa$, $\lambda$, $\mu$. There is a nice formula for $\kappa$: $$ \begin{align*} \kappa \;&=\; \int_{-\infty}^\infty |g_3(x)|\,dx \\[6pt] &=\; e^{-\alpha^2/2}(\alpha^2-\alpha-3)\,-\,e^{-\beta^2/2}(\beta^2-\beta-3)\,+\,e^{-\gamma^2/2}(\gamma^2-\gamma-3) \\[6pt] &\approx\; 3.5519079 \end{align*} $$ where $\alpha<\beta<\gamma$ are the three roots of the polynomial $x^3 - x^2 - 5x + 1$.

The formula for $\lambda$ is similarly nice: $$ \begin{align*} \lambda \;&=\; \int_{-\infty}^\infty \frac{g_3(x) g_4(x)}{|g_3(x)|}dx \\[6pt] &=\; e^{-\alpha^2/2}p(\alpha)\,-\,e^{-\beta^2/2}p(\beta)\,+\,e^{-\gamma^2/2}p(\gamma) \\[6pt] &\approx\; -3.307248 \end{align*} $$ where $p(x) = \dfrac{1}{12}\bigl(3x^4-34x^2+12x+25\bigr)$.

Things get a little bit dicey after that, since the values of $\psi_1'(0)$, $\psi_2'(0)$, and $\psi_3'(0)$ come into play. In particular, $$ \begin{align*} \mu \;&=\; \int_{-\infty}^\infty \frac{g_3(x) g_5(x)}{|g_3(x)|}dx \,+\, \bigl(2g_4(\alpha)+g_3'(\alpha)\bigr)\psi_1'(0) \\[6pt] &\qquad -\, \bigl(2g_4(\beta)+g_3'(\beta)\bigr)\psi_2'(0) +\, \bigl(2g_4(\gamma)+g_3'(\gamma)\bigr)\psi_3'(0) \end{align*} $$ It's possible to compute the values of $\psi_1'(0)$, $\psi_2'(0)$, and $\psi_3'(0)$ by examining the gradient of $f(x,\epsilon)/\epsilon^3$ near the points $(\alpha,0)$, $(\beta,0)$, and $\gamma(0)$. The result is that $$ \psi_1'(0) = q(\alpha),\qquad \psi_2'(0)=q(\beta),\qquad \psi_3'(0)=q(\gamma) $$ where $$ q(x) \;=\; \frac{3x^5-46x^3+12x^2+93x-12}{12(x^4-x^3-8x^2+3x+5)}. $$ In particular, $$ \psi_1'(0)\approx -0.832825,\qquad \psi_2'(0) \approx 0.0971987,\qquad \psi_3'(0)\approx 1.06896.$$ Using these formulas, I'm getting that $\mu\approx 3.70537$, but this is complicated enough that I'm not very confident about this value.

Edit: In the comments, Clement asks how we know that the remainder term in the expansion $$ \Delta(\epsilon) \;=\; \kappa \epsilon^3 \,+\, \lambda \epsilon^4 \,+\,\mu\epsilon^5 \,+\, o(\epsilon^5) $$ is indeed $o(\epsilon^5)$. Well, observe that the function $f(x,\epsilon)$ has a simple antiderivative with respect to $x$: $$ \begin{align*} F(x,\epsilon) \;&=\; \epsilon\sqrt{\frac{\pi}{2}}\left(\mathrm{erf}\left(\frac{x-\epsilon}{\sqrt2}\right)-\mathrm{erf}\left(\frac{x}{\sqrt2}\right)\right) \\ &\qquad+\, (1-\epsilon)\sqrt{\frac{\pi}{2}}\left(\mathrm{erf}\left(\frac{x}{\sqrt{2(1+\epsilon)}}\right)-\mathrm{erf}\left(\frac{x-\epsilon^2}{\sqrt{2(1+\epsilon)}}\right)\right) \end{align*} $$ where $\mathrm{erf}$ is the error function. Since $\mathrm{erf}$ is an entire function, it is clear that $F(x,\epsilon)$ is analytic on $\mathbb{R}\times(-1,1)$. It is easy to check that $F(x,\epsilon)\to 0$ as $x\to\pm\infty$, so $$ \Delta(\epsilon) \;=\; 2 F(\psi_1(\epsilon),\epsilon)-2F(\psi_2(\epsilon),\epsilon)+2F(\psi_3(\epsilon),\epsilon) $$ Since $\psi_1$, $\psi_2$, and $\psi_3$ are $C^\infty$, we can compute the power series for $\Delta(\epsilon)$ in the usual way, giving the results above. The remainder is $o(\epsilon^5)$ because of Taylor's Theorem.

(Edited. I had the wrong sign between the two main parts of $f$.)

The idea is: Get a global overview and then develop anything in sight into a power series with respect to $t$ (which is your $\epsilon$).

Begin by replacing $f$ with the function $$g(x,t):=e^{x^2/2} f(x,t)\ .$$ Then $g$ can be written in the neighborhood of $(0,0)$ in the form $$g(x,t)={1\over2}(-1+5x+x^2-x^3) t^3 +{1\over24}(12 -93 x+\ldots) t^4+\ldots\ .$$ The polynomial $-1+5x+x^2-x^3$ has three real zeros $x_1\doteq-1.9$, $x_2\doteq0.195$, $x_3\doteq2.7$. These will be the constant terms in three functions $t\mapsto x=\psi_i(t)$ that describe where $x\mapsto f(x,t)$ changes sign for fixed small $t$.