The spherical symmetric solution in your question provides a basis for the PDE in the infinite space $R^n$. Your boundary value problem in a sphere of a finite radius, you need to expand the solution in the series of products of the spherical Bessel function of $r$ and $n$-dimensional spherical harmonics.
The Laplace operator in $3$-dimensional spherical coordinate is
$$0=-\frac{\partial u}{\partial t}+\Delta u =\frac{\partial u}{\partial t}-\frac{1}{r^2}\frac{\partial}{\partial r} \left(r^2 \frac{\partial u}{\partial r}\right) + \frac{1}{r^2 \sin\theta} \frac{\partial}{\partial \theta} \left(\sin\theta \frac{\partial u}{\partial \theta}\right) + \frac{1}{r^2 \sin^2\theta} \frac{\partial^2 u}{\partial \varphi^2} \tag1$$
Suppose $u(r,\theta,\phi)=T(t)R(r)\Theta(\theta)\Phi(\varphi)$. Substitute it into Equation (1) and divide it by $f$, we have
$$0=-\frac1T\frac{\partial T}{\partial t}+\frac{1}{r^2R}\frac{\partial}{\partial r} \left(r^2 \frac{\partial R}{\partial r}\right) + \frac{1}{r^2\Theta \sin\theta} \frac{\partial}{\partial \theta} \left(\sin\theta \frac{\partial \Theta}{\partial \theta}\right) + \frac{1}{r^2 \sin^2\theta\Phi} \frac{\partial^2 \Phi}{\partial \varphi^2} \tag2$$
We see that the first term is independent of $(r,\theta,\varphi)$ while the remaining terms are independent of $t$. The first term must then be a constant. Let us call it $-\lambda$. We have
$$T'(t)=-\lambda T(t) \tag3$$
Multiplying (2) by $r^2$, we see that the third and fourth terms are independent of $r$ while the second term depends only on $r$. So we introduce another constant $-\mu$.
$$\frac{1}{r^2R}\frac{d}{d r} \left(r^2 \frac{d R}{d r}\right)+\Big(\lambda-\frac\mu {r^2}\Big)R=0. \tag4$$
Apply the same trick again, we obtain
$$\Phi''(\varphi)+\nu\Phi(\varphi)=0 \tag5$$
for some constant $\nu$ and for $\Theta(\theta)$
$$\Theta''+\cot\theta\,\Phi'+(\mu-\nu(\csc\theta)^2)\Theta=0 \tag6$$
for some constant $\mu$.
We will solve Equation ($4$) for with boundary condition $R(r=1)=0$ where we set the radius of the sphere at $1$.
...
Jumping to the conclusion first, let
$$q_{\lambda,k,m}(r,\theta,\varphi):=j_k(r\sqrt\lambda)Y_{k,m}(\theta,\varphi),\ \lambda>0,$$
where $j_k$ is the spherical Bessel function with $\sqrt\lambda$ being a positive root of it, $Y_{k,m}$ is the spherical harmonic function. Expand
$$g(r,\theta,\varphi)=\sum_{\lambda,|m|\le k}a_{\lambda,k,m}q_{\lambda,k,m}(r,\theta,\varphi)$$
to solve for the coefficient $a_{\lambda,k,m}$. This is facilitated by the property that the set $\{q_{\lambda,k,m}(r,\theta,\varphi)\}_{\lambda,k,m}$ is orthogonal in $L^2$.
In the case of the homogeneous PDE where $f=0$,
$$u(r,\theta,\varphi,t) = \sum_{\lambda,|m|\le k}a_{\lambda,k,m}q_{\lambda,k,m}(r,\theta,\varphi)e^{-\lambda t}.$$
By Duhamel's principle, the inhomogeneous equation can be viewed as linear superposition of the homogeneous problem and
$$u(r,\theta,\varphi,t) = \sum_{\lambda,|m|\le k}a_{\lambda,k,m}q_{\lambda,k,m}(r,\theta,\varphi,t)e^{-\lambda t}+\int_0^t ds \sum_{\lambda,|m|\le k}a_{\lambda,k,m}(s)q_{\lambda,k,m}(r,\theta,\varphi)e^{-\lambda (t-s)},$$
where
$$f(r,\theta,\varphi,t)=\sum_{\lambda,|m|\le k}a_{\lambda,k,m}(t)q_{\lambda,k,m}(r,\theta,\varphi).$$
This is how one builds up the solution in a sphere. More details can be found in e.g. ON GREEN'S FUNCTIONS IN THE THEORY OF HEAT
CONDUCTION IN SPHERICAL COORDINATES by Arnold N. Lowan and Green’s Function for the Heat Equation by Abdelgabar Adam Hassan.
Now if $g(x)$ is spherically symmetric thus is only a radial function $g(r)$, then the solution is simpler.
rn−1w′+12rnw=A.
Evans goes on to state that as r→∞, w,w′→0 so that A=0. This can be seen from the form of the solution to the ODE. If A≠0 then w,w′ cannot tend to zero as r grows.
– asd Sep 23 '18 at 13:39