I'm aware this is an old question and may be considered somewhat obsolete, but let me answer it nevertheless - if only for the sake of completeness.
The position space representation of the Klein-Gordon Green function (propagator) clearly looks intimidating. The trick is to do the calculation in momentum space, where the propagator is just a rational function. Before actually doing this let me point out that in the massless case, $m = 0$, and for a static source, $\rho = \rho (\mathbf{x})$, one is actually solving the Poisson equation. If the source is radially symmetric, $\rho = \rho(r)$, as suggested in the question, the solution is the Coulomb potential, $\phi = \phi(r) \sim 1/r$. Allowing for non-vanishing mass, one obtains a Yukawa potential, $\phi \sim \exp(-mr)/r$.
This can be shown explicitly in terms of Fourier integrals. First, transform the field,
$$ \phi(k) = \int d^4 x \, e^{i k\cdot x} \, \phi (x) \; , $$
and similarly the the source, $\rho \to \rho(k)$. The momentum space solution of the KG equation is then
$$ \phi(k) = \frac{\rho(k)}{k^2 - m^2 + i \epsilon} \; , $$
with $k^2 = k_0^2 - \mathbf{k}^2$ and the causal $i\epsilon$-prescription. (Change appropriately for retarded/advanced solutions.) Then transform back to position space,
$$ \phi (x) = \int \frac{d^4 k}{(2\pi)^4} \, e^{-i k\cdot x} \, \frac{\rho(k)}{k^2 - m^2 + i \epsilon} \; . \quad (**)$$
As an example, take a Gaussian source, $\rho(r) = \rho_0 \exp(-\alpha r^2)$, with 'normalisation' $ \rho_0 = (\alpha/\pi)^{3/2}$. Its Fourier transfrom is $\rho(k) = 2\pi i \, \delta(k^0) \exp (-k^2/4\alpha)$. The $k^0$-integral in $(**)$ is hence trivial, and the $d^3 k$ integration can be done with the usual residue technique writing $\mathbf{k}^2 + m^2 \equiv (\kappa+im)(\kappa-im)$. The result is
$$ \phi(r) = e^{m^2/4\alpha} \, \frac{e^{-mr}}{4\pi r} \; .$$
In the point source limit ($\alpha \to \infty$) we reobtain the standard Yukawa potential.
For time dependent sources there will be energy transfer (no $\delta(k^0)$), and the integral $(**)$ will normally be harder. Typically, one can do the $k^0$-integration via residues and the remaining one(s) using stationary phase as e.g. in Ch. 2.1 of Peskin and Schroeder.
Green's functions are not unique. Any solution of that satisfies the homogeneous equation, $$(\partial_t^2 - \nabla^2 + m^2)f = 0$$ in the region of interest can be added to the Green's function without spoiling the inhomogeneous equation. The homogeneous part must, therefore, be set by boundary conditions. You can get the Bessel function part of the Feynman version of the Green's function by solving for the stationary Green's function with a 4-dimensional space, that is:
$$\left(-\frac{1}{r^3} \frac{\partial }{\partial r} \left[r^3 \frac{\partial}{\partial r}\right] + m^2\right)G(r) = \frac{i}{\Omega_4 r^3}\delta(r)$$ and analytically continuing to imaginary time. Inside of the forward light cone, where the argument of $K$ is imaginary, gives a Hankel function of the first kind: $$\begin{align}\frac{m}{4\pi^2\sqrt{\tau^2}}K_1\left(im\sqrt{\tau^2}\right) &= -\frac{m}{8\pi\sqrt{\tau^2}} H_1^{(2)}\left(m\sqrt{\tau^2}\right) \\
& = -\frac{m}{8\pi\sqrt{\tau^2}} \left[J_1\left(m\sqrt{\tau^2}\right) - i Y_1\left(m\sqrt{\tau^2}\right)\right],\end{align}$$ with $J_a$ and $Y_a$ the Bessel functions of the first and second kind, respectively.
Because both the operator and the inhomogeneous part are real, the imaginary part of the Green's function must be a solution to the homogeneous equation and the real part must solve the inhomogeneous. That is:
$$\left(\frac{\partial}{\partial t}^2 - \nabla^2 + m^2\right) \operatorname{Im}\left\{\frac{m}{4\pi^2\sqrt{\tau^2}}K_1\left(im\sqrt{\tau^2}\right)\right\} = 0,\ \mathrm{and} \\
\left(\frac{\partial}{\partial t}^2 - \nabla^2 + m^2\right) \operatorname{Re}\left\{\frac{m}{4\pi^2\sqrt{\tau^2}}K_1\left(im\sqrt{\tau^2}\right)\right\} = -\delta^4(x),$$ so the Green's function is the real part of that expression.
Features that remain to be shown at this point in this post: that the split into retarded and advanced propagators is valid (needed to preserve causality), and the zero mass limit of the propagator gives the correct limit. The energy injected by the impulse cannot be tracked because it is infinite.
Best Answer
You have found the reason that most 4D QFTs need to have renormalization: The propagator is (UV) divergent!
To "cure" the theory, you need to regularize the theory (make the divergences expressible in some simple parameter like a momentum cutoff) - e.g. by dimensional regularization - and then proceed with some renormalization scheme. That the theory is non-interacting means that there are only finitely many divergent things - the propagator itself - as opposed to the infinitely many divergent diagrams what would need to be organized in orders of perturbation theory in an interacting theory, but, as you see, it does not preclude the necessity of renormalization as such.