If $f$ is well behaved, the solution $u(t)$ can be obtained analytically. First, define $w(t,x) \equiv u(t,x)+F(x)$ where $\Delta F(x)=f$. Assuming $f$ is not pathological, such $F$ exists and can be calculated directly, because we know the kernel of Laplace's equation:
$F(x)=\int_{\mathbb{R}^n} \frac{f(y)}{|x-y|}dy $
This $w$ obeys a homogeneous equation:
$\partial_t w= \partial_t (u+F)=\partial_t u=\Delta u+f=\Delta (u+F)= \Delta w$
The solution for $w$ can also be obtained because we know the kernel of the heat equation:
$w(t,x)=\frac{1}{(4\pi t)^{n/2}} \int_{\mathbb{R}^n} w(0,y)e^{-\frac{|x-y|^2}{4t}}dy$
It is seen from this solution that if $w(x,0)$ decays fast enough with x, then $w\to 0$ for $t\to\infty$ and for all $x$. Therefore, $\Delta u \to -f$ as needed.