By the Helffer-Sjöstrand formula,
$$
f(-\Delta_D) = \frac1\pi \int_{\mathbb C} (\partial_{\bar z} \tilde
f)(z) (-\Delta_D -z)^{-1}\,L(dz),
$$
where $L(dz)$ is the Lebesgue measure and $\tilde f$ is an almost-analytic extension of $f$. Now the resolvent $(-\Delta_D-z)^{-1}$ for $z\in\mathbb C\setminus[0,\infty)$ has kernel
$$
K_z(x,y) = \begin{cases}
\frac1{2\mu}\left(e^{-\mu(x-y)}-e^{-\mu(x+y)}\right), & 0< y< x, \\
\frac1{2\mu} \left(e^{\mu(x-y)}-e^{-\mu(x+y)}\right), & x< y< \infty,
\end{cases}
$$
where $\mu=\sqrt{-z}$ with $\Re\mu>0$. After a short computation this reveals
\begin{align*}
K_f(x,y) &= \frac1\pi \,\int_{\mathbb C} (\partial_{\bar z}\tilde f)(z)
K_z(x,y)\,L(dz) \\ &= \frac1\pi \int_0^\infty
f(\lambda)\,\frac{\sin(\sqrt{\lambda} x)\sin(\sqrt{\lambda}
y)}{\sqrt{\lambda}}\,d\lambda
\end{align*}
for the kernel of $f(-\Delta_D)$. (An alternative approach uses the spectral resolution of $-\Delta_D$, see e.g. on page 261 here.)
If you want $f(-\Delta_D)\colon H^{-n}(\mathbb R_+) \to L^2(\mathbb R_+)$ for $n\in \mathbb N$, then by duality you need $\bar f(-\Delta_D)\colon L^2(\mathbb R_+) \to H_0^n(\mathbb R_+)$. However, $$ (\partial_x K_{\bar f})(0,y) = \frac1\pi \int_0^\infty \bar f(\lambda)\sin(\sqrt{\lambda} y)\,d\lambda, $$ which is not identically zero (except when $f\equiv0$$f|_{\,\mathbb R_+}\equiv0$). So, you already fail to map $L^2(\mathbb R_+)$ into $H_0^2(\mathbb R_+)$.