Basic integration over Dirac delta For positive $x$ and continuous support function $f(y,z)$, is it correct to conclude that,
$$\lim_{\varepsilon\rightarrow0^{+}}\int_{0}^{x}dy\int_{0}^{y}dz\,f(y,z)\,\delta(y-\varepsilon)\,\delta(z-\varepsilon)\,=\,f(0,0)$$
Is there any way to simplify the result for general value of $x$ and $\varepsilon$?
 A: We will see that the problem is about the value of the Dirac delta at the boundary of an integral.
When the Dirac is not on the boundary there is no problems, that is when $c\neq a$ and $c\neq b$ you could write
$$\label{1}\tag{1}
\int_a^b g(x)\,\delta_c(\mathrm d x) = g(c)\,\mathbf 1_{(a,b)}(c)
$$
where I used the notation $\delta_a(\mathrm  d x)$ instead of $\delta(x-a)\,\mathrm d x$.
Now, mathematically speaking, the Dirac delta is a measure for which discrete sets can have positive measure. Hence one should differentiate
$$\label{2}\tag{2}
\int_{[a,b]} f(x)\, \delta_a(\mathrm d x) = f(a)
$$
and
$$\label{3}\tag{3}
\int_{(a,b)} f(x)\, \delta_a(\mathrm d x) = 0
$$
Hence in these cases, the notation $\int_a^b$ is not clear. Now in problems coming from physics, the Dirac delta could come from some limit, in this case one should come back to the limit, and not simplify it to quickly as a Dirac delta. For instance, when the Dirac delta was a limit of a very narrow positive even function $g_n$ of integral $1$, such as a Gaussian, then one will rather obtain
$$\tag{4}\label{4}
\int_a^b f(x)\, g_n(x-a)\,\mathrm d x \underset{n\to\infty}{\to} f(a)/2.
$$
Hence, the value of \eqref{1} when $c=a$ or $c=b$ might depend on your problem.
Now, what is the link with your problem? Your integral can be written
$$
I_\varepsilon = \int_0^x\int_0^y f(y,z) \,\delta_{\varepsilon}(\mathrm d z) \,\delta_{\varepsilon}(\mathrm d y)
$$
Starting with the inner integral and taking $\varepsilon > 0$, it holds
$$
\int_0^y f(y,z) \,\delta_{\varepsilon}(\mathrm d z) = f(y,\varepsilon) \,\mathbf 1_{\varepsilon<y}(\varepsilon),
$$
without clear meaning when $\varepsilon =y$.
Hence, the full expression for which you want to take the limit is
$$
I_\varepsilon = \int_0^x f(y,\varepsilon) \mathbf 1_{\varepsilon < y}\,\delta_{\varepsilon}(\mathrm d y) = f(\varepsilon,\varepsilon) \,\mathbf 1_{\varepsilon <\varepsilon},
$$
so we are exactly in the case where we are missing the information of what we should do about the boundary value: what I wrote as $1_{\varepsilon <\varepsilon}$ is just the unknown value that would be $1$ in case \eqref{2}, would be $0$ in case \eqref{3} and would be $1/2$ in case \eqref{4}. Other values would be possible if your Dirac delta was coming from a specific type of approximating functions.
If $f$ is continuous and the boundary value of your Dirac delta was giving some value $c(\varepsilon)$ depending continuously of $\varepsilon$, then the limit will be
$$
c(0)\, f(0,0).
$$
