Engineering PDE with dynamic BC and source term The PDE:
$$\frac1D C_t-Q=\frac2rC_r+C_{rr}$$
on the domain $r \in [0,\bar{R}]$ and $t \in [0,+\infty]$ and where $D$ and $Q$ are Real constants. We're looking for a function $C(r,t)$.
The BC:
$$C(0,t)=f(t)$$
$$C_r(\bar{R},t)=0$$
The IC:
$$C(r,0)=C_0$$
If $f(t)=0$ then I know the solution. Assume:
$$C(r,t)=C_E(r)+v(r,t)$$
$$-Q=\frac2rC_r+C_{rr}$$
$$-Q=\frac2rC_E'(r)+C_E''(r)$$
$$rC_E''+2C_E'+Qr=0$$
where $C_E(r)$ is the steady-state solution ($t \to \infty$).
$$C_E(r)=-\frac{Qr^2}{6}+\frac{c_1}{r}+c_2$$
$$\text{because }\lim_{r\to 0}C(r)=\infty\Rightarrow c_1=0$$
But since as $f(t) \neq 0$, $c_2$ cannot be determined.
All help will be appreciated.

Edit. In the case of $f(t)=0$ the solution, summarised, becomes:
$$c_2=0$$
$$C_E(r)=-\frac{Qr^2}{6}$$
$$C(r,t)=-\frac{Qr^2}{6}+v(r,t)$$
Compute partial derivatives:
$$C_t=v_t$$
$$C_r=-\frac{Qr}{3}+v_r$$
$$C_{rr}=-\frac{Q}{3}+v_{rr}$$
Inserting in the PDE then gives the homogeneous PDE in $v(r,t)$:
$$\frac1D v_t=\frac2r v_r+v_{rr}$$
Ansatz: $v(r,t)=R(r)T(t)$, then separation of variables yields the ODE solutions, with $-m^2$ a separation constant:
$$T(t)=c_3\exp(-m^2 D t)$$
$$R(r)=c_4\frac{\sin mr}{r}$$
BCs:
$$R(0)=0$$
$$R'(\bar{R})=0$$
$$R'=c_4\frac{mr\cos mr-\sin mr}{r^2}$$
$$R'(\bar{R})=c_4\frac{m\bar{R}\cos m\bar{R}-\sin m\bar{R}}{\bar{R}^2}=0$$
The eigenvalues $m_i$ are the solutions to the transcendental equation:
$$m_i\bar{R}=\tan m_i\bar{R}$$
So we have:
$$v(r,t)=\sum_{i=1}^\infty A_i\exp(-m_i^2 D t)\frac{\sin m_ir}{r}$$
Determine the $A_i$ the usual way with the IC and the Fourier series.
So we have:
$$C(r,t)=-\frac{Qr^2}{6}+\sum_{i=1}^\infty A_i\exp(-m_i^2 D t)\frac{\sin m_ir}{r}$$
 A: NB: I am using the $(t,r)$ convention as used by Olver, etc.
You seem to be working with the inhomogeneous heat equation, which is a well known problem. Via a change of spacial and temporal variables we can reduce your equation to
$$(\partial_t-\Delta)u=q$$
Where $q$ is constant. We couple this with the initial and boundary conditions
$$u(0,r)=u_0\\ u(t,0)=f(t)\\ (\partial_r u)(t,R)=0$$
Where again $u_0$ is constant. The usual method one would use to do this is by Fourier transform, but since we are working in the domain $r\in[0,R]$ and not $r\in[0,\infty)$, it's not very fruitful. So let's come up with a different approach. Let's say we've found a solution to the IBVP
$$(\partial_t-\Delta)u=0\\u(0,r)=u_0\\u(t,0)=f(t)\\(\partial_ru)(t,R)=0\tag{1}$$
Call it $u_\mathrm{h}$. And let's now suppose we've found a solution to the IBVP
$$(\partial_t-\Delta)u=q\\u(0,r)=0\\u(t,0)=0\\(\partial_r u)(t,R)=0$$
Call it $u_\mathrm{p}$ (it seems you've already found a formula for it as well, just set $C_0=0$). Now if we consider the function
$$v=u_{\mathrm h}+u_{\mathrm p}$$
Since the operator $\partial_t-\Delta$ is linear, we can see that
$$(\partial_t-\Delta)v=q\\v(0,r)=u_0\\u(t,0)=f(t)\\(\partial_rv)(t,R)=0$$
In other words, $v$ is the solution of your problem.
So all you need to do now is solve $(1)$. Since the RHS is zero, separation of variables should work nicely here. There might be concerns with uniqueness, etc, here, but since this problem arises in engineering, I doubt you'll be too worried.
A: Solving $\boldsymbol{(1)}$ from my previous answer
I've adapted my solution method from this link. The trick is to take a transform in the temporal variable. We can use either a one sided Fourier transform or a Laplace transform, it doesn't really matter. I will go with the Laplace transform as it avoids complex numbers. Recall that our problem was
$$(\partial_t-\Delta)u=0\\ u(0,r)=u_0\\u(t,0)=f(t)\\(\partial_r u)(t,1)=0$$
Laplace transform in time:
$$\mathcal{L}_t(\partial_tu)-\mathcal{L}_t(\Delta u)=0$$
We use the fact that the Laplace transform and the spacial derivatives commute, and the simple identity
$$\mathcal{L}_t(\partial_t u)(s,r)=s~(\mathcal{L}_tu)(s,r)-u(0,r)$$
(To prove this, just treat $r$ as a constant and use integration by parts.)
Using now the abbreviation $\mathcal{L}_tu=U$ and $u(0,r)=u_0$ we get the equation
$$(\Delta U)(s,r)-s~U(s,r)=-u_0$$
Treating $s$ as constant, this is just the Helmholtz equation for $U$, with an extra term on the right hand side. We can deal with this extra term using a particular solution. We can see that
$$\frac{u_0}{s}$$
Is a particular solution of the above, and combining this with my approach in my other answer we add this to the homogeneous solution, getting
$$U(s,r)=\frac{u_0}{s}+A(s)\frac{e^{-\sqrt{s}~r}}{r}+B(s)\frac{e^{\sqrt{s}~r}}{r}$$
Using the third condition, (and of course because the Laplace transform of $0$ is $0$) we get
$$B(s)=A(s)e^{-2\sqrt{s}}\frac{\sqrt{s}+1}{\sqrt{s}-1}$$
So then
$$U(s,r)=\frac{u_0}{s}+A(s)\left(\frac{e^{-\sqrt{s}~r}}{r}+e^{-2\sqrt{s}}\frac{\sqrt{s}+1}{\sqrt{s}-1}\frac{e^{\sqrt{s}~r}}{r}\right)$$
Using the fact that the ILT of $u_0/s$ is just $u_0$, and the convolution theorem, we get
$$u(t,r)=u_0+\int_0^\infty a(\tau)\phi(t-\tau,r)\mathrm d\tau$$
Where
$$ a:=\mathcal{L}^{-1}(A)\\ \phi:=\mathcal{L}_s^{-1}\left((s,r)\mapsto \frac{e^{-\sqrt{s}~r}}{r}+e^{-2\sqrt{s}}\frac{\sqrt{s}+1}{\sqrt{s}-1}\frac{e^{\sqrt{s}~r}}{r}\right)$$
In practice $\phi$ can be calculated numerically. Unfortunately, the initial and boundary data don't transform nicely. We get
$$u(0,r)=u_0\implies U(0,r)=\mathcal{L}(u_0)(0)=\text{undefined}$$
Since the Laplace transform of a constant, say $c$, is $c/s$. One can also see that trying to evaluate $U(s,0)$ doesn't work either. So unfortunately, to work out $a$, one needs to solve the integral equation
$$u(t,0)=f(t)=u_0+\int_0^\infty a(\tau)\phi(t-\tau,0)\mathrm d\tau$$
For $a$. $\phi$ is in principle already known, so we can solve for $a$ numerically, and then we know $u(t,r)$ in its entirety.
