$\renewcommand{\Re}{\operatorname{Re}}$ $\newcommand{\erf}{\operatorname{erf}}$ $\newcommand{\dag}{\dagger}$ $\newcommand{\const}{\mathrm{const}}$
Let $n\ge 3$. Consider spherically symmetric density $f$ (thus depending on $r=(x_1^2+x_2^2+\ldots+x_n^2)^{\frac{1}{2}}$ only). Then it creates a density which is also spherically symmetric.
Assume first that $f=0$ in $B(0,R)$ and consider $u$ in $B(0,R)$. Here $u$ must be harmonic and then due to mean-value theorem $u=u(0)$ is constant (in particular $\nabla u=0$). More presisely \begin{equation} u= -\frac{1}{n-2} \int_R^\infty \rho f(\rho)\,d\rho \label{eq-1} \end{equation} where we replaced lower limit $0$ by $R$ since $f(\rho)=0$ as $r< R$.
Assume now that $f=0$ in $\{r\ge R\}$ and consider $u$ there. Then $u=r^{2-n} A+B$ (Home Assignment 8, problem 1) and to have potential at infinity equal to $0$ we must take $B=0$. Further plug it into (26.2) with $\Omega=B(0,R)$. The left-hand expression becomes the total charge $\int_{B(0,R)} f(x)\,dV$ while the right-hand expression becomes $-(n-2)\sigma_n A$ and therefore \begin{equation} u= -\frac{1}{(n-2)\sigma_n}r^{2-n}\int _{B(0,R)} f\,dx = -\frac{1}{n-2}r^{2-n}\int _0^R f(\rho)\,d\rho. \label{eq-2} \end{equation} Then in the general case we to calculate $u$ we break $f=f_1+f_2$ where $f_1=0$ for $r\ge R$ and $f_2=0$ for $r\le R$ and then calculate $u$ as $r=R$ and finally set $R=r$: \begin{equation} u= -\frac{1}{n-2} \int_r^\infty \rho f(\rho)\,d\rho-\frac{1}{n-2}r^{2-n}\int _0^r f(\rho)\,d\rho. \label{eq-3} \end{equation}
\begin{multline} u(y)=\int_{\Omega} G^0(x,y)\Delta u(x)\,dV + \\ \int_\Sigma \bigl(-u(x)\frac{\partial G^0}{\partial \nu_x}(x,y)+G^0(x,y)\frac{\partial u}{\partial \nu}(x)\bigr)\, dS \qquad \tag{26.7}\label{eq-26-7} \end{multline} with \begin{equation} G^0(x,y)=\left\{\begin{aligned} -&\frac{1}{(n-2)\sigma_n}|x-y|^{2-n}&& n\ne 2,\\ -&\frac{1}{4\pi}|x-y|^{-1} &&n=3,\\ &\frac{1}{2\pi}\log |x-y| &&n=2,\\ &\frac{1}{2} |x-y| &&n=1 \end{aligned}\right. \tag{26.8} \end{equation} where we changed notation $G^0$ instead of $G$ as we will redefine $G$ later.
Let $g(x,y)$ be solution to problem
\begin{align}
&\Delta_x g(x,y)=0\qquad&&\text{in }\ \Omega,\label{eq-4}\\[3pt]
&g(x,y)=-G^0(x,y) &&\text{as }\ x\in\Sigma\label{eq-5}
\end{align}
and condition $u\to 0$ as $|x|\to \infty$ if $\Omega$ is unbounded. In virtue of (\ref{eq-4}) and
(26.5)
\begin{multline}
0=\int_{\Omega} g(x,y)\Delta u(x)\,dV + \\
\int_\Sigma \bigl(-u(x)\frac{\partial g}{\partial \nu_x}(x,y)+g(x,y)\frac{\partial u}{\partial \nu}(x)\bigr)\, dS. \qquad
\label{eq-6}
\end{multline}
Adding to (\ref{eq-26-7}) we get
\begin{multline}
u(y)=\int_{\Omega} G(x,y)\Delta u(x)\,dV + \\
\int_\Sigma \bigl(-u(x)\frac{\partial G}{\partial \nu_x}(x,y)+G(x,y)\frac{\partial u}{\partial \nu}(x)\bigr)\, dS \qquad
\label{eq-7}
\end{multline}
with
\begin{equation}
G(x,y):=G^0(x,y)+g(x,y).
\label{eq-8}
\end{equation}
So far we have not used (\ref{eq-5}) which is equivalent to $G(x,y)=0$ as $x\in \Sigma$. But then
\begin{equation*}
u(y)=\int_{\Omega} G(x,y)\Delta u(x)\,dV
- \int_\Sigma \frac{\partial G}{\partial \nu_x}(x,y)u(x)\, dS
\end{equation*}
and therefore
\begin{equation}
u(y)=\int_{\Omega} G(x,y)f(x)\,dV
- \int_\Sigma \frac{\partial G}{\partial \nu_x}(x,y)\phi(x)\, dS
\label{eq-9}
\end{equation}
is a solution to
\begin{align}
&\Delta u=f \qquad&&\text{in }\ \Omega,\label{eq-10}\\[3pt]
&u=\phi &&\text{on }\ \Sigma.\label{eq-11}
\end{align}
Similarly if we replace (\ref{eq-5}) by Robin boundary condition \begin{equation} \bigl(\frac{\partial g}{\partial \nu_x}-\alpha g\bigr)(x,y)=-\bigl(\frac{\partial G^0}{\partial \nu_x}-\alpha G^0\bigr) (x,y) \qquad \text{as }\ x\in\Sigma \label{eq-12} \end{equation} we get $\bigl(\frac{\partial G}{\partial \nu_x}-\alpha G\bigr)(x,y)=0$ as $x\in \Sigma$ and rewriting (\ref{eq-7}) as \begin{multline} u(y)=\int_{\Omega} G(x,y)\Delta u(x)\,dV + \\ \int_\Sigma \bigl(-u(x) \bigl[ \frac{\partial G}{\partial \nu_x}-\alpha G \bigr] (x,y)+G(x,y)\bigl[ \frac{\partial u}{\partial \nu}-\alpha u\bigr]\bigr)\, dS \qquad \label{eq-13} \end{multline} we get \begin{equation} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_\Sigma G(x,y)\psi(x)\, dS \label{eq-14} \end{equation} for solution of the (\ref{eq-10}) with the boundary condition \begin{equation} \frac{\partial u}{\partial \nu}-\alpha u=\psi \qquad \text{on }\ \Sigma.\label{eq-15} \end{equation}
Finally, if we take $g$ satisfying mixed boundary condition \begin{align} & g(x,y)=-G^0(x,y) && \text{as }\ x\in \Sigma',\label{eq-16}\\ & \bigl(\frac{\partial g}{\partial \nu_x}-\alpha g\bigr)(x,y)=-\bigl(\frac{\partial G^0}{\partial \nu_x}-\alpha G^0\bigr) (x,y) &&\text{as }\ x\in\Sigma''\qquad\label{eq-17} \end{align} we get \begin{multline} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_{\Sigma'} G(x,y)\psi(x)\, dS - \\ \int_{\Sigma''} \frac{\partial G(x,y)}{\partial \nu_x}\phi(x)\, dS\qquad \label{eq-18} \end{multline} for solution of the (\ref{eq-10}) with the boundary condition \begin{align} &u=\phi && \text{on }\ \Sigma',\label{eq-19}\\ &\frac{\partial u}{\partial \nu}-\alpha u=\psi && \text{on }\ \Sigma'',\label{eq-20} \end{align} where $\Sigma=\Sigma'\cup\Sigma''$ and $\Sigma'\cap\Sigma''=\emptyset$.
Definition 1. $G(x,y)$ (defined for corresponding boundary problem) is called Green's function.
Remark 1. It is similar (by usage) to heat kernel $\frac{1}{2\sqrt{kt}}e^{-\frac{|x-y|^2}{4kt}}$.
One can prove
Theorem 1. \begin{equation}G(x,y)=G(y,x). \label{eq-21} \end{equation}
Consider now purely Neumann problem in the connected domain. We cannot solve $\Delta g=0$ with the boundary condition $\frac{\partial g}{\partial \nu_x}=-\frac{\partial G^0}{\partial \nu_x}$ as such problem requires one solvability condition to the right-hand expression and boundary value \begin{equation} \int_\Omega f\, dV+ \int _\Sigma \psi\,dS =0 \label{eq-22} \end{equation} and this condition fails as $\int_\Sigma \frac{\partial G^0}{\partial \nu_x}\, dS\ne 0$ (one can prove it).
To fix it we consider $g$ satisfying \begin{align} &\Delta_x g(x,y)=c\qquad&&\text{in }\ \Omega,\label{eq-23}\\[3pt] &\frac{\partial g}{\partial \nu_x}(x,y)=-\frac{\partial G^0}{\partial \nu_x}(x,y) &&\text{as }\ x\in\Sigma\label{eq-24} \end{align} with unknown constant $c$; chosing it correctly one can satisfy solvability condition. Then \begin{equation*} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_{\Sigma} G(x,y)\psi(x)\, dS + c \int _V u\,dx \end{equation*} and therefore \begin{equation} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_{\Sigma} G(x,y)\psi(x)\, dS +C \label{eq-25} \end{equation} gives us solution if it exists (and it exists under solvability condition (\ref{eq-24})). This solution is defined up to a constant.