The proof is essentially the same, the only difference is that with $\varphi \in \mathscr{S}(\mathbb{R}^3)$ we can't use one fixed outer radius $R$ of the spherical shell over which we integrate, since Schwartz functions generally don't have compact support. But Schwartz functions decay rapidly as $\lVert x\rVert \to \infty$, and thus if we let the radius of the outer sphere in Green's formula tend to $\infty$, that part of the boundary integral tends to $0$.
Since $\varphi \in \mathscr{S}(\mathbb{R}^3) \implies \nabla^2 \varphi \in \mathscr{S}(\mathbb{R}^3)$ and $f(x) = \lVert x-x_0\rVert^{-1}$ is a locally integrable function of polynomial growth - in fact, $\lim\limits_{\lVert x\rVert\to\infty} f(x) = 0$ - we have
$$\int_{\mathbb{R}^3} f(x) \nabla^2\varphi(x)\,d\mu = \lim_{\substack{\varepsilon \to 0 \\ R \to \infty}} \int_{A(\varepsilon,R;x_0)} f(x)\nabla^2\varphi(x)\,d\mu,\tag{1}$$
where $A(\varepsilon,R;x_0) = \{ x\in \mathbb{R}^3 : \varepsilon < \lVert x-x_0\rVert < R\}$ and $0 < \varepsilon < R < \infty$. In $(1)$, we can take iterated limits in either order, or a simultaneous limit. Since $f\cdot \nabla^2\varphi$ is integrable, all these limits are well-defined and coincide.
Now we use $\nabla^2 f(x) = 0$ on $\mathbb{R}^3 \setminus \{x_0\}$ and Green's formula to rewrite
\begin{align}
\int_{A(\varepsilon,R;x_0)} f(x)\nabla^2\varphi(x)\,d\mu
&= \int_{A(\varepsilon, R; x_0)} f(x)\nabla^2\varphi(x) - \varphi(x)\nabla^2 f(x)\,d\mu \\
&= \int_{\partial A(\varepsilon, R; x_0)} f(x) \frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial \nu}(x)\,dS,
\end{align}
where $\frac{\partial}{\partial \nu}$ is the normal derivative (in direction of the outer normal).
$\partial A(\varepsilon, R; x_0)$ consists of two pieces, the outer sphere $S_R = \{ x : \lVert x-x_0\rVert = R\}$ and the inner sphere $S_{\varepsilon} = \{ x : \lVert x-x_0\rVert = \varepsilon\}$. For the integral over $S_{\varepsilon}$, the argument that
$$\lim_{\varepsilon \to 0} \int_{S_{\varepsilon}} f(x) \frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial \nu}(x)\,dS = -4\pi \varphi(x_0)$$
is completely identical to the case of compactly supported $\varphi$. It remains to see that
$$\lim_{R \to \infty} \int_{S_R} f(x) \frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial \nu}(x)\,dS = 0.\tag{2}$$
But that follows from the fast decay of $\varphi$ and its derivatives. By definition of $\mathscr{S}(\mathbb{R}^3)$, for all $N \in \mathbb{N}$ there is a constant $C \in (0,+\infty)$ such that $\lvert \lVert x\rVert^N\varphi(x)\rvert \leqslant C$ and $\bigl\lvert \lVert x\rVert^N\frac{\partial \varphi}{\partial x_k}(x)\bigr\rvert \leqslant C$ for all $x\in \mathbb{R}^3$ and $1\leqslant k \leqslant 3$. For $R > 2\lVert x_0$ we have
$$\lVert x\rVert = \lVert (x-x_0) + x_0\rVert \geqslant \lVert x-x_0\rVert - \lVert x_0\rVert = R - \lVert x_0\rVert \geqslant \frac{R}{2}$$
on $S_R$, and so we can estimate
$$\lvert\varphi(x)\rvert \leqslant \frac{C}{\lVert x\rVert^N} \leqslant \frac{2^NC}{R^N},\quad \biggl\lvert \frac{\partial \varphi}{\partial\nu}(x)\biggr\rvert \leqslant \frac{2^N\sqrt{3}\,C}{R^N},$$
which yields
$$\Biggl\lvert \int_{S_R} f(x) \frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial \nu}(x)\,dS \Biggr\rvert \leqslant \biggl(\frac{1}{R}\cdot \frac{2^N\sqrt{3}\,C}{R^N} + \frac{2^NC}{R^N}\cdot \frac{1}{R^2}\biggr)\cdot 4\pi R^2 \leqslant \frac{\tilde{C}}{R^{N-1}},$$
and we see that taking any $N > 1$ gives us $(2)$.
It may be worth noting that the argument isn't specific to dimension $3$. For any dimension $d \neq 2$, the function $f_d(x) = \lVert x-x_0\rVert^{2-d}$ is harmonic on $\mathbb{R}^d\setminus \{x_0\}$, and the same computation shows $\nabla^2 f_d = (2 - d)\omega_{d-1}\cdot \delta_{x_0}$, where $\omega_{d-1}$ is the $d-1$-dimensional volume of the unit sphere in $\mathbb{R}^d$. The last estimate becomes
$$\Biggl\lvert \int_{S_R} f(x) \frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial \nu}(x)\,dS \Biggr\rvert \leqslant \biggl(\frac{1}{R^{d-2}}\cdot \frac{2^N\sqrt{d}\,C}{R^N} + \frac{2^NC}{R^N}\cdot \frac{\lvert d-2\rvert}{R^{d-1}}\biggr)\cdot \omega_{d-1} R^{d-1} \leqslant \frac{\tilde{C}}{R^{N-1}}.$$
For dimension $d = 2$, the corresponding function is $f_2(x) = -\log \lVert x-x_0\rVert$ with $\nabla^2 f_2 = 2\pi \delta_{x_0}$ and for the integral over the outer boundary circle we eventually get an estimate $\tilde{C}\cdot \frac{\log R}{R^{N-1}}$.
