You have an $\epsilon$ on the right, but no $\epsilon$ on the left, so that's weird. Anyway, suppose $f(x,y)=\frac{1}{4\pi}\ln(x^2+y^2)$ for $(x,y)\neq(0,0)$ (define it however you wish at the origin). Then, $f$ is Locally Lebesuge-integrable on $\Bbb{R}^2$ so it defines a distribution. For each $\epsilon>0$, let $f^{\epsilon}(x,y)=\frac{1}{4\pi}\ln(\epsilon^2+x^2+y^2)$. Then, $f_{\epsilon}$ is smooth. By the dominated convergence theorem, $f^{\epsilon}\to f$ in the topology of distributions $\mathcal{D}'(\Bbb{R}^2)$. Therefore, it also follows that all distributional derivatives converge; i.e for any multindex $\partial^{\alpha}(f^{\epsilon})\to \partial^{\alpha}f$ in the topology of $\mathcal{D}'(\Bbb{R}^2)$. In particular, the Laplacian $\Delta f^{\epsilon}\to \Delta f$ in the topology of $\mathcal{D}'(\Bbb{R}^2)$
So, we now calculate the Laplacian of $f^{\epsilon}$ (which because $f^{\epsilon}$ is smooth, its standard Laplacian coincides with the distributional one):
\begin{align}
(\Delta f^{\epsilon})(x,y)&=\frac{1}{\pi}\frac{\epsilon^2}{(\epsilon^2+x^2+y^2)^2}
\end{align}
Now, define $\zeta(x,y)=\frac{1}{\pi}\frac{1}{(1+x^2+y^2)^2}$ and $\zeta_{\epsilon}(x,y):=\frac{1}{\epsilon^2}\zeta\left(\frac{x}{\epsilon},\frac{y}{\epsilon}\right)$, Then, the above calculation shows that $\Delta f^{\epsilon}=\zeta_{\epsilon}$. Since $\zeta\in L^1(\Bbb{R}^2)$, and in fact $\int_{\Bbb{R}^2}\zeta=1$ (use polar coordinates to verify this), it follows that $\zeta_{\epsilon}\to \delta$ in $\mathcal{D}'$ as $\epsilon\to 0^+$. Thus, putting everything together,
\begin{align}
\Delta f&=\lim_{\epsilon\to 0^+}\Delta (f^{\epsilon})=\lim_{\epsilon\to 0^+}\zeta_{\epsilon}=\delta.
\end{align}
