Solving the diffusion equation to model an injection of a drug: IC $c(r,t=0)=c_0\ 0\leq r \leq r_0$ I'm trying to model an injection of a substance into a fluid, assuming no flow is created. I'm considering the following:
$$\frac{\partial c(r,t)}{\partial t} = 
\frac{1}{r^2} 
\frac{\partial }{\partial r}
\left(
r^2
\frac{\partial }{\partial r}
c(r,t)
\right),
\quad r\in\mathbb{R^3},t>0$$
$$c(r,0)=c_0\quad 0\leq r \leq r_0$$
$$c(r,0)=0\quad r>0$$
Where $c=c(r,t)$ is the concentration of the substance being injected. Initially $c=c_0>0$ in a sphere centred at the origin, the injection is therefore assumed to occur instantly and for the substance to form a perfect sphere around the point of injection.
My question is how do we solve this system giving us the behavior of the sphere as time moves forward. I'm most interested in very shortly after $t=0$, therefore I've been thinking of this as being solved in $\mathbb{R^3}$, but if it simplifies the problem $r$ can lie within an arbitrarily sized sphere with radius $R_0$ centred at the origin, with no flux conditions. 
It's been awhile since I've worked with PDEs so would appreciate any advice or references to where this problem has been solved before.
 A: The answer depends a lot on how you want to implement your model: either in all of $\mathbb{R}^3$ or in a finite ball of radius $R>0$.  
Case 1 - All of $\mathbb{R}^3$
In this case you want to solve
$$
\begin{cases}
\partial_t u(x,t) = \Delta u(x,t) & \text{ for } x \in \mathbb{R}^3, t >0 \\
u(x,0) = u_0(x),
\end{cases}
$$ 
where 
$$
u_0(x) = 
\begin{cases}
c_0 & \text{if }|x| \le r_0 \\
0 &\text{otherwise.}
\end{cases}
$$
This can be solved with the usual heat kernel via convolution:
$$
u(x,t) = \int_{\mathbb{R}^3} \frac{1}{(4\pi t)^{3/2}} e^{-|x-y|^2/(4t)} u_0(y)dy = \frac{c_0}{(4\pi t)^{3/2}} \int_{B(0,r_0)} e^{-|x-y|^2/(4t)} dy.
$$
The resulting solution will be radial, bounded, smooth, and decay at infinity.
Case 2 - The ball $B(0,R)$
Obviously we need $R>r_0$ here.  In this case the key thing is that you need to choose boundary conditions to enforce on $\partial B(0,R)$, i.e. for $x \in \mathbb{R}^3 $ such that $|x|=R$.  The two standard choices would be the Dirichlet condition, where you specify the value of $u$, and the Neumann condition, where you specify the value of $\nabla u \cdot \nu$ (the normal derivative).  Based on what you're trying to model, I would guess that the Neumann condition is more appropriate.  If you specify that the solution satisfies $\nabla u \cdot \nu =0$ then this means that there is no "flux of the drug" out of the domain, which is quite reasonable.  Moreover, since your initial data is radial this can be written as $\partial_r u =0$ for $|x|= R$.
You thus want to solve the $3D$ radial heat equation:
$$
\begin{cases}
\partial_t u(r,t) = \partial_r^2 u(r,t) + \frac{2}{r} \partial_r u(r,t) & \text{for }0 < r < R, t>0 \\
\partial_r u(R,t) = 0 & \text{for } t >0 \\
u(r,0) = u_0.
\end{cases}
$$
This problem can be solved by using separation of variables.  You assume an ansatz $u(r,t) = S(r) T(t)$ and plug in and use ODE theory to find an infinite family of solutions.  Then  you represent the solution to your problem as an appropriate linear combination of these solutions.  For details I will refer you to this post.  In it they consider different initial conditions, but the analysis of the separation of variables is the same.  The post also restricts to $R=1$, so you will need to modify things a bit to handle general $R$.  
