I have very little background with functional derivatives and I would like to clarify some issues.

I am trying to compute the second functional derivative of the Klein Gordon action expressed in real components

$S[\phi_1,\phi_2,\partial_{\mu}\phi_1,\partial_{\nu}\phi_2]=\int{}d^4x\,\mathcal{L}=\int{}d^4x(-\partial_{\mu}\phi_1\partial^{\mu}\phi_1-\partial_{\mu}\phi_2\partial^{\mu}\phi_2+$

$-m^2(\phi_1^2+\phi_2^2)-\lambda(\phi_1^2+\phi_2^2)^2)$

I know how to compute the first functional derivative

$$\frac{\delta{}S}{\delta\phi_1(x)}=\frac{\partial\mathcal{L}}{\partial\phi_1}-\partial_{\mu}\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\phi_1)}$$

which gives

$$2\left(\partial^2\phi_1(x)-m^2\phi_1(x)-4\lambda\phi_1(x)[\phi_1^2(x)+\phi_2^2(x)]\right)$$

from now on I have doubts. I wanna compute

$$\frac{\delta^2S}{\delta\phi_1(y)\delta\phi_1(x)}$$

If there were no derivatives this would be trivial but I don't know how to deal with $\partial^2$. I have thought in introducing a delta to write this as an integral and use the formula I have used to compute the first functional derivative, but I am concerned about taking derivatives of something having a Dirac delta.

So, my question is, how am I supposed to compute this functional derivative?