In THIS ANSWER and THIS ONE, I provided primers on the Dirac Delta.
The notation $\int_a^b f(x)\delta(x-c)\,dx$ is interpreted to mean the functional $\langle fp_{ab},\delta_c\rangle$.
Here, $p_{ab}$ is the "rectangular pulse" function, $p_{ab}(x)=u(x-a)-u(x-b)$, and $u$ is the unit step (or Heaviside Function) where
$$u(x)=\begin{cases}1&,x>0\\\\0&,x<0\end{cases}$$
Note that there are various conventions for the value $u(0)$.
Therefore, with $a=\epsilon$, $b=-\epsilon$, and $c=0$, we have
$$\begin{align}
\int_{-\epsilon}^\epsilon f(x)\delta(x)\,dx&=\langle fp_{-\epsilon\,\epsilon},\delta_0\rangle\\\\
&=\int_{-\infty}^\infty p_{-\epsilon\,\epsilon}(x)f(x)\delta(x)\,dx \\\\
&=f(0)p_{-\epsilon\,\epsilon}(0)\\\\
&=f(0)
\end{align}$$
as was to be shown!