This is a similar question to [this](https://mathoverflow.net/questions/470070/estimating-a-solution-to-an-euler-type-ode/470111?noredirect=1#comment1220446_470111) but with a different ODE. 

Let $f$ be a continuous function in $L^2([1,\infty)$ satisfying $\sup_{r\geq 1} r|f(r)| <\infty$. Let $\ell$ be a positive integer, $R>1$, and $a$ be a real number. 

Let $u(r)$ be a function on $[R,\infty)$ solving the following IVP: 
$$(r^2-1)u''(r)+2ru'(r) - \ell(\ell+1) u(r) = f(r), \quad r\in [R,\infty)$$
$$u(R)= a, \quad \lim_{r\to \infty} u(r) = 0$$

We should be able to estimate the solution $u$ by $f$ and $a$. I want to understand precisely how $\ell$ will appear in this estimate.

I am speculating the following: there exists a constant $C$ independent of $\ell$, $a$, $u$ and $f$ satisfying 
$$\ell \sup_{r\geq R}r|u(r)| + \sup_{r\geq R}r^2|u'(r)| \leq C(\ell^{-1}\sup_{r\geq R} r|f(r)| + \ell^{-1}\lVert f\rVert_{L^2} + \ell|a|)$$

Can anyone verify if this is true? How would I prove this estimate (or the correct one in case the above is incorrect)? 

The solutions to the homogeneous problem are the Legendre polynomials of the first and second kind $P_{\ell}$, $Q_{\ell}$. We can use the method variation of parameters to get

$$u(r) = AQ_{\ell}(r) + P_{\ell}(r) \int_{r}^{\infty} Q_{\ell}(t) f(t)W(t)^{-1}(t^2-1)^{-1}dt + Q_{\ell}(r) \int_{R}^r P_{\ell}(t) f(t) W(t)^{-1}(t^2-1)^{-1} dt$$

where $W = P_{\ell}Q_{\ell}'-P_{\ell}'Q_{\ell}$ and $A$ is some constant that will depend on $a$. I have tried to use this formula to get the estimate but couldn't. 

Any help is appreciated. 


EDIT: Using the Frobenius method we get,  

$$P_{\ell}(r) = \sum_{n=0}^{\ell} a_n r^{\ell-n}, \quad Q_{\ell}(r) = \sum_{n=0}^{\infty}b_n r^{-\ell-1-n}$$
where $a_n$, $b_n$ are recursively defined by $a_0 = b_0 = 1$ (by choice), $a_1 = b_1 = 0$, and for $n\geq 2$
$$a_n = \frac{(\ell-n+1)(\ell-n+2)}{n^2-n(2\ell+1)} a_{n-2}, \quad b_n = \frac{(\ell+n-1)(\ell+n)}{n(n+2\ell+1)}b_{n-2}$$

Then the Wronskian satisfies: 
$$W(t)(t^2-1) = 2\ell+1$$
in which I used Lagrange's identity to deduce that $W(t)(t^2-1)$ is a constant and took the limit as $t$ goes to $\infty$ to find the constant. 

We would then like to get the following estimates: 

$$\left| rP_{\ell}(r) \int_{r}^{\infty} Q_{\ell}(t) f(t)dt \right|\leq \frac{C}{\ell} \sup_{s\geq R}s|f(s)|$$
$$\left| rQ_{\ell}(r) \int_{R}^r P_{\ell}(t) f(t)dt \right|\leq \frac{C}{\ell} \sup_{s\geq R}s|f(s)|$$
for some $C$ independent of $\ell$. This will then imply $ \sup_{r\geq 1}r|u(r)| \leq C\ell^{-2}\sup_{r\geq R} r|f(r)|$, which is part of what we're trying to show.(Note that the second factor of $\ell$ in the right side comes from $W$).

It will be sufficient if we have estimates on $Q_{\ell}$ and $P_{\ell}$ looking like: 

$$(1) \qquad \sup_{r\geq R}r^{\ell+1}|Q_{\ell}(r)| \leq C, \quad \sup_{r\geq R}r^{-\ell} |P_{\ell}(r)| \leq C$$ 

Where $C$ is independent of $\ell$ (it can depend on $R$). If so, then the estimate I'm trying to get follows using the same calculation @IgorKhavkine performed. 

However, I haven't been able to get those estimates in equation (1). 

I think I could use the recursive relation: 
$$(\ell+1)P_{\ell+1}(r) = (2n+1)r P_{\ell}(r) - \ell P_{\ell-1}$$
to get the desired estimate on $P_{\ell}$. 

I have tried to make use of the identity: 
$$Q_{\ell}(r) = B\int_{0}^{\infty} (r+(r^2-1)^{1/2}\cosh \theta)^{-\ell-1} d\theta$$
where $B^{-1} = \int_{0}^{\infty} (1+\cosh \theta)^{-\ell-1} d\theta$ ( added so that it agrees with the expression of $Q_{\ell}$ written earlier). But I couldn't make use of this to get the desired estimate on $Q_{\ell}$. I would need to bound the expression 
$$\frac{\int_{0}^{\infty} (1+(1-R^{-2})^{1/2}\cosh \theta)^{-\ell-1} d\theta}{\int_{0}^{\infty} (1+\cosh \theta)^{-\ell-1} d\theta}$$
by a constant independent of $\ell$. But it seems too difficult.