What conditions must satisfy the source term of an ODE given some boundary conditions? Suppose you have this second order linear ordinary differential equation:
$$ \frac{\mathrm d}{\mathrm d x}\left[\left(1-\frac{2}{x}\right)\frac{\mathrm du(x)}{\mathrm d x}\right]-(l+1)(l+2) \frac{u(x)}{x^2}=S(x)$$
where $u(x)$ is the unknown, $S(x)$ is a source term and $l$ is an integer $l\geq0$. (It is a particular case of the Jacobi equation)
This equation has two regular singular points in $x=0,2,\infty$. I am looking for a solution defined in $x\in [2,\infty)$ so that $u(x)$ is regular in $x=2$ and:
$$\lim_{x\to\infty}u(x)=0$$
$S(x)$ is finite in $x=2$ and goes to $0$ for $x\to\infty$. 
To satisfy the second boundary condition I found (using Mathematica) that you need to have for $S(x)$ the asymptotic behaviour$^1$:
$$S(x)\sim \frac{1}{x^3},\qquad(\text{for }x\to\infty) $$
Why is this the case?
More in general, given a non-homogeneous ODE and Dirichlet-like boundary conditions (like in my example), what should the behaviour of $S$ be on the boundary points so that the equation has a solution?

NOTES
$1)$ Here's how I found this:


*

*Set $S(x) = \frac{1}{x}$ and $l$ to some integer number 

*Solved the equation with Mathematica getting two integration constants $c_1$ and $c_2$  

*Expanded in series the solution around $\infty$  

*The solution is a polynomial in $x$ (+ $\mathcal O(x^{-1})$ terms) where the coefficient of the $x$ term and the constant term are independent of those constants. So no matter what you can't get a convergent solution for $r\to\infty$

*Repeat for $S(x) = \frac{1}{x^2}$ to get a solution that is constant at infinity and for $S(x) = \frac{1}{x^3}$ to get a solution that goes to $0$ as wanted.

 A: The general question for arbitrary ODE is quite difficult. But in this specific case, there is some estimate of $S(x)$ that you can obtain.
Note that $\lim_{x\rightarrow \infty} u(x) = 0$, thus for large $x$, your equation is approximated by
$$x^2 \frac{d}{dx}\left[\left(1 - \frac 2x\right)\frac{du(x)}{dx}\right] = S(x)$$
Which is equivalent to
$$\frac{d}{dx}\left[\left(1 - \frac 2x\right)\frac{du(x)}{dx}\right] = \frac{S(x)}{x^2}$$
Since $\lim_{x\rightarrow \infty} u(x) = 0$, the limit of its derivative must behave similarly, so $\lim_{x\rightarrow \infty} u'(x) = 0$, by barbalat's lemma. Consequently, the LHS is $0$. Thus,
$$\lim_{x\rightarrow\infty} \frac{S(x)}{x^2} = 0$$
This is true only if $S(x)$ has order smaller than $x^2$, e.g. $S(x) = x$, etc. I can't see precisely how Mathematica comes up with $S(x) \sim \frac{1}{x^3}$ for large $x$, but that does satisfy the above condition.

A very informal analysis does gives a better asymptotic estimation though. Since $\lim_{x\rightarrow\infty} u(x) = 0$, suppose $u(x)$ is of order $\frac{1}{x}$. This means for large $x$, $u'(x)$ is of order $\frac{1}{x^2}$. "Patching" things together you have something like this
$$x^2 \frac{d}{dx}\left[O\left(\frac{1}{x}\right) O\left(\frac{1}{x^2}\right)\right] - O\left(\frac{1}{x}\right) = S(x)$$
This is, loosely speaking, equivalent to
$$O\left(x^2\right) O\left(\frac{1}{x^4}\right) - O\left(\frac{1}{x}\right) = S(x)$$
Or
$$O\left(\frac{1}{x^2}\right) - O\left(\frac{1}{x}\right) = S(x)$$
Since $S(x)$ drives the equation, it has to goes to $0$ faster than any term on the left hand side, thus $S(x)$ is of order $\frac{1}{x^3}$ or lower. If it is of a lower order than $\frac{1}{x^3}$, then we can expect $u(x)$ to be of a lower order. But since $u(x) \sim \frac{1}{x}$ is perfectly fine assumption (without solving the equation and this is a very informal analysis), we conclude that
$$S(x) \sim \frac{1}{x^3}$$

Edit: with the new information, then I think my very informal analysis makes more sense (since you were simply trying out different solution for $S(x)$).
