In the equation
$$R''(r)+\frac5rR'(r)+\frac4{r^2}LR(r)+2ER(r)-\frac2rUR(r)=0\tag1$$
$L$ and $U$ are constant matrices, and $E$ is a scalar. I'm looking for a solution, which decays to $0$ as $r\to\infty$. I've been able to find asymptotic solution looking like
$$R(r)\sim\exp\left(-\sqrt{-2E}r\right)\;\text{as }r\to\infty.\tag2$$
But I'd like to have at least one $U$-dependent term so as to know ratios of the components of $R$ at infinity. I tried to multiply by $r$ both sides of the equation and take $r\to\infty$, but this didn't work because the $R''$ and $2ER$ terms appeared to be dominant and went to infinity (these are the terms which helped me get $(2)$).
So how can I get the next term in the expansion of the solution at infinity?
We have a nice fact about the terms with $L$ and $U$: one of the terms is much smaller than the other, and moreover, it's $L$-dependent term which goes to zero faster:
$$\frac4{r^2}L-\frac2rU\sim-\frac2rU\text{ as }r\to\infty.\tag{I}$$
This gives us the relation with $U$ but without $L$:
$$R''(r)+\frac5rR'(r)-\frac2rUR(r)\sim-2ER(r)\text{ as }r\to\infty.\tag{II}$$
Now we have only one matrix, $U$, and this allows us to switch to a basis where it's diagonal. If we denote
$$U=KWK^{-1},\tag{III}$$
where $W$ is a diagonal matrix with $u_i$ on the diagonal ($i=1,2,...,n$ for $n\times n$ matrix $U$), then, denoting
$$Q(r)=K^{-1}R(r),\tag{IV}$$
we can transform the relation $(\mathrm{II})$ to
$$Q''(r)+\frac5rQ'(r)-\frac2rWQ(r)\sim-2EQ(r)\text{ as }r\to\infty.\tag{V}$$
But since $W$ is diagonal, we can actually split this system into a set of independent relations for components $Q_i$ of $Q$:
$$Q_i''(r)+\frac5rQ_i'(r)-\frac2r u_i Q_i(r)\sim-2EQ_i(r)\text{ as }r\to\infty.\tag{VI}$$
These relations can be treated by the usual method of dominant balance as follows. Set
$$Q_i(r)=\exp(S_i(r)).\tag{VII}$$
We'll get
$$(S_i'(r))^2+S_i''(r)+\frac5r S_i'(r)\sim\frac2ru_i-2E.\tag{VIII}$$
Assuming $(S_i'(r))^2\gg S_i''(r)$ as $r\to\infty$, we drop $S_i''$ term, simplifying $(\mathrm{VIII})$ to
$$(S_i'(r))^2+\frac5r S_i'(r)\sim\frac2ru_i-2E.\tag{IX}$$
Solving quadratic equation for $S_i'(r)$ and integrating, we get
$$S_i(r)\sim-\frac52\ln r\pm\int\limits_1^r\sqrt{\frac{25}{4r^2}+\frac2ru_i-2E}\,\mathrm dr+C_i.\tag X$$
As we want the decaying solution, we choose $\pm\to-$. Thus we have a set of $Q_i$, each with its own arbitrary constant factor of $e^{C_i}$. Going back to original basis as given by $(\mathrm{IV})$, we finally get our $R(r)$, up to $n$ arbitrary constants, which should be defined from other considerations.