Consider the following partial differential equation $$ \frac{\partial f (r,t)}{\partial t} - \frac{k}{r^2} \left( \frac{\partial}{\partial r} \left( r^ \, \frac{\partial f (r,t)}{\partial r} \right) - 2f(r,t) \right) + \frac{f(r=1, t)}{r^2} = 0\, , $$ subject to the boundary conditions $$ \left. \frac{\partial f}{\partial r} \right|_{r=1} = 0\, , \qquad f(r=\infty, t) = 0 \, , $$ together with the initial condition $f(r, t=0) = \epsilon$. Here, $k$ and $\epsilon$ are positive real numbers, with $|\epsilon| \ll 1$
To search an analytical solution of this PDE, i used the method of separation of variables. Accordingly, the solution is expressed as $f(r,t) = F(r)G(t)$. Using this approach, it can be shown that $G(t)$ satisfies the following differential equation $$ G'(t) + \operatorname{const.} G(t) = 0 \, . $$ The latter admit the exponential function as a solution. The problem is that, by considering an decaying exponential, one would obtain a vanishing solution in the long-time limit as $t \to \infty$. This is in contradiction with the numerical solution where, for instance, $f(r=1,t)$ reached a steady/plateau value in the limit $t \to \infty$. Does this mean that one cannot employ the method of separation of variables in this problem?
I'm assuming that the equation in the title is the correct one and that in the text has some typo.
The general solution to $$\frac{\partial f (r,t)}{\partial t} - \frac{k}{r^2} \left( \frac{\partial}{\partial r} \left( r^ \, \frac{\partial f (r,t)}{\partial r} \right) - 2f(r,t) \right) + g(r) = 0\, ,$$ is of the form $$f(r,t)=H(r,t)+P(r,t),$$ where $P(r,t)$ is any particular solution of the equation and $H(r,t)$ is the general solution of the homogeneous equation:$$\frac{\partial f (r,t)}{\partial t} - \frac{k}{r^2} \left( \frac{\partial}{\partial r} \left( r^ \, \frac{\partial f (r,t)}{\partial r} \right) - 2f(r,t) \right) = 0\ .$$
My guess is that you are forgoting the particular solution $P(r,t)$. A particular $P$ can be find assuming that it's a function of $r$ only. With this ansatz we get: $$ - \frac{k}{r^2} \left( \frac{d}{dr} \left( r^ \, \frac{dP}{dr} \right) - 2P \right) + g(r) = 0$$ wich simplifies to $$rP''+P'-2P=\frac{r^2g(r)}{k},$$ whoose solutions can be found with modified Bessel functions and using the variation of parameters. Depending on your's $g(r)$ the form of the $P$ can be really ugly. But there's no contradiction. As you find, the $H(r,t)$ is a sum of terms $F(r)G(t)$ where the G's are decaying exponentials, so in the long-time limit your full solution $H+P$ tends to $P(r)$.