SO the heat equation is
$$\partial_t u=\lambda\nabla^2u$$
With $\lambda$ being some diffusivity constant. We can use the Laplacian in spherical coordinates and remove the $\theta,\phi$ dependency to get
$$\partial_t u=\frac{\lambda}{r^2}\partial_r\left(r^2\partial_ru\right)$$
Let's assume the equation is separable.
$$\frac{G'(t)}{G(t)}=\frac{\lambda}{r^2F(r)}\partial_r\left(r^2F'(r)\right)=-c$$
Let's focus on the spacial part. Make the substitution $F(r)=f(r)/r$ so
$$F'(r)=\frac{r\cdot f'(r)-f(r)}{r^2}$$
Which leads us to
$$f''(r)=\frac{-c}{\lambda}f(r)$$
$$f(r)=a_1\sin\left(\sqrt{\frac{c}{\lambda}}r\right)+a_2\cos\left(\sqrt{\frac{c}{\lambda}}r\right)$$
Imposing the boundary conditions, namely $f(0)=f(a)=0$ we require $a_2=0$ and we get constants
$$c_n=\lambda \frac{n^2\pi^2}{a^2}$$
Which gives our spacial eigenfunctions
$$F_n(r)=A_n\operatorname{sinc}\left(\frac{n\pi}{a}r\right)=A_nj_0\left(\frac{n\pi}{a}r\right)=A_n\sqrt{\frac{a}{2nr}}J_{1/2}\left(\frac{n\pi}{a}r\right)$$
Using Bessel functions (this will be important later)
The time component follows:
$$\frac{G'(t)}{G(t)}=-c=-\lambda \frac{n^2\pi^2}{a^2}$$
$$G_n(t)=B_n\exp\left(-\lambda \frac{n^2\pi^2}{a^2}t\right)$$
So we can write our solution as (using $A_n\sqrt{\frac{a}{2n}}\cdot B_n= C_n$)
$$u(r,t)=\sum_{n=1}^\infty C_n\frac{J_{1/2}\left(\frac{n\pi}{a}r\right)}{\sqrt{r}}\exp\left(-\lambda \frac{n^2\pi^2}{a^2}t\right)$$
Initial conditions:
We can see that
$$\sqrt{r}\cdot u(r,0)=\sum_{n=1}^\infty C_nJ_{1/2}\left(\frac{n\pi}{a}r\right)$$
This is a perfect time to use Fourier-Bessel series. Let $\alpha_{1/2,n}$ be the $n$th root of $J_{1/2}$, which in this case is the rather simple $\alpha_{1/2,n}=n\pi$.
So
$$\sqrt{r}u_0(r)=\sum_{n=1}^\infty C_nJ_{1/2}\left(\frac{\alpha_{1/2,n}}{a}r\right)$$
The literature then tells us that
$$C_n=\frac{\int_0^a r^{3/2}u_0(r)J_{1/2}\left(\frac{\alpha_{1/2,n}}{a}r\right)\mathrm{d}r}{(a/2)J_{3/2}(\alpha_{1/2,n})^2}$$
EXAMPLE: Let $u_0(r)=(r^2-1)\log(r)$. The first few coefficients are
$$\{C_n\}=\{0.624023,0.386661,0.10893,...\}$$
Here is a plot of the initial data $u_0(r)$ as well as its first five partial series-

And the Mathematica code used (so you can try yourself):
Subscript[u, 0][r_] := (r^2 - 1) Log[r]
c[n_] := 2*
NIntegrate[
r^(3/2) Subscript[u, 0][r] BesselJ[1/2, n*\[Pi]*r], {r, 0, 1}]/
BesselJ[3/2, n*\[Pi]]^2
Subscript[u, 0][r_, N_] :=
Sum[c[n]*(Sqrt[2] Sin[n*\[Pi]*r])/(\[Pi]*r*Sqrt[n]), {n, 1, N}]
Plot[{Subscript[u, 0][r], Subscript[u, 0][r, 1],
Subscript[u, 0][r, 2], Subscript[u, 0][r, 3], Subscript[u, 0][r, 4],
Subscript[u, 0][r, 5]}, {r, 0.01, 1},
PlotLegends -> {"\!\(\*SubscriptBox[\(u\), \(0\)]\)(r)",
"\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;1)",
"\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;2)",
"\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;3)",
"\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;4)",
"\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;5)"}]
Here $u_0(r;N)$ is the $N$th partial series. Another example when $u_0(r)=\sqrt{1-r}$ :

EDIT: What if we have Neumann, not Dirichlet, boundary conditions? That is, what if we want
$$\partial_r[u](0,t)=\partial_r[u](a,t)=0$$
Well, let's return to our spacial part of our solution:
$$F(r)=a_1\frac{\sin\left(\sqrt{\frac{c}{\lambda}}r\right)}{r}+a_2\frac{\cos\left(\sqrt{\frac{c}{\lambda}}r\right)}{r}$$
The derivative is
$$F'(r)=\frac{\cos\left(\sqrt{\frac{c}{\lambda}}r\right)\left(a_1\sqrt{\frac{c}{\lambda}}r-a_2\right)-\sin\left(\sqrt{\frac{c}{\lambda}}r\right)\left(a_1+a_2\sqrt{\frac{c}{\lambda}}r\right)}{r^2}$$
In order to get $F'(0)=0$ we require $a_2=0$. Then our expression for the derivative simplifies to
$$F'(r)=a_1\frac{\sqrt{\frac{c}{\lambda}} r\cos\left(\sqrt{\frac{c}{\lambda}}r\right)-\sin\left(\sqrt{\frac{c}{\lambda}}r\right)}{r^2}$$
What we're trying to do is find values of $c$ such that $F'(a)=0$ given $\lambda,a$. Unfortunately, the roots of $F(r)$ are not easily obtainable in closed form. Letting $\beta_{\lambda,n}$ be the $n$th positive root of $F'(r)$ given the constant $\lambda$ (we can find these via numerical algorithms) we require
$$\beta_n=\sqrt{\frac{c_n}{\lambda}}a\implies c_n=\lambda \left(\frac{\beta_n}{a}\right)^2$$
So our spacial eigenfunctions are
$$F_n(r)=A_n\operatorname{sinc}\left(\frac{\beta_n}{a}r\right)=A_n j_0\left(\frac{\beta_n}{a}r\right)=A_n\sqrt{\frac{\pi}{2\frac{\beta_n}{a}r}}J_{1/2}\left(\frac{\beta_n}{a}r\right)\to A_n\frac{1}{\sqrt{r}}J_{1/2}\left(\frac{\beta_n}{a}r\right)$$
Where of course we redefine $A_n$ as we go along for convenience. After this point, we follow the same process as before with the Fourier Bessel series and whatnot - just in this case we have to compute the roots $\beta_n$ beforehand.