your calculation is fine but there is no need to employ numerical methods for the given equation at hand. The equation can be solved using separation of variables with a physical ansatz for the radial part. Of course you need not restrict your problem to a radial symmetry of the solution, the full Laplace will lead to spherical harmonics for the angular parts. The following derivation can directly be applied in this case.
Solution using separation of variables
Taking $$T(r,t) = f(t)\cdot B(r)$$
results in
$$f'(t) B(r) = \alpha f(t)\frac{1}{r^2}\partial_r \left(r^2 B'(r) \right)$$
Now we see that we can re-write this as
$$\frac{f'(t)}{f(t)} = \alpha \frac{1}{r^2}\frac{\partial_r \left(r^2 B'(r) \right)}{B(r)}$$
and see that both sides are independent of each other because of the different dependencies, hence they must be a constant, say $-\lambda$, for the left hand side $$\frac{f'(t)}{f(t)} = -\lambda$$ for which we find the solution
$$f(t) = f_0 e^{-\lambda t}\ .$$
$\lambda$ can be seen as a dissipation constant defining the loss of heat per time unit. I found a nice illustration of this process on youtube or on wikipedia for another geometry:
Now to the radial part,
$$\alpha \frac{1}{r^2}\partial_r \left(r^2 B'(r) \right) = -\lambda B(r)$$
for which we find
$$B(r) = c_1 \frac{e^{-\mathrm{i}\kappa r}}{r} + c_2\frac{e^{\mathrm{i}\kappa r}}{r} $$
with $$\kappa = \sqrt{\lambda/ \alpha}\ .$$
Q: Can you tell me what kind of solutions we just arrived at?
Application to a spherical body
Now this was the general solution to the problem but we have to apply it to your situation. First of all, wee need some initial thermal distribution: $$T_0(r) = T(r,0)$$ to determine the dynamics from this point. Second of all, we need to know what it means to have some sphere in our system. We don't have a choice but to assign some distribution to $\alpha$, determining different heat coefficients in different regions, e.g.
$$\alpha = \begin{cases}
\alpha_{i} & r\leq R\\
\alpha_{0} & r>R\end{cases} $$
The only thing left to do is to determine the constants $\lambda$ and $c_i$ from the continuity of $T$ at $r=R$ and meaningful values at $r = \infty$ ($T\rightarrow 0$) and $r = 0$ ($T$ finite).
I will leave this step to you.
Having the solution at hand you can compare it to your simulation results and further extend it, e.g. to a system of shells or special boundary conditions.
Sincerely
