Asymptotic behavior of nonlinear boundary value problem ODE solutions I am trying to solve numerically the following real-valued system $(f,h)$ on $[0,\infty]$:
$$-h''-\frac{1}{r}h'+\lambda_c h(f^2-1)+4\lambda_h h^3=0$$
$$-f''-\frac{1}{r}f'+\frac{1}{r^2}f+ f(-\lambda_f+2\lambda_c h^2)+\lambda_f f^3=0$$
where $\lambda_h, \lambda_f, \lambda_c $ are positive constants, and with boundary conditions $f(0)=0, f(\infty)=1, h'(0)=0, h(\infty)=0$.  The condition $h'(0)$ is necessary for $h$ to be finite at the origin. We may assume without loss of generality $h(0)>0$, and for simplicity we can set $\lambda_f=\lambda_h=1$.
In order to solve this system, it is helpful to know the asymptotic behavior of $(f,g)$ at $r=0$ and $r=\infty$.
My question is: what is the asymptotic behavior of $(f,h)$ as $r\to \infty$.
In particular, if an asymptotic expansion at infinity exists for $(f,h)$, what are e.g the first two terms ?

I find that at small $r$, $f(r) \sim \alpha r $ and $h(r) \sim a +br^2 +cr^4$ for constants $\alpha, a,b,c$.
At large $r$, if there is no coupling between $f$ and $h$ (i.e $\lambda_c=0$) then it's easy to show $f(r)=1-C / r^2+\mathcal{O}\left(1 / r^4\right)$.
Assuming $f$ has the same form when the coupling is turned on, and assuming $h$ admits a  large-$r$ expansion in powers of $1/r$, I find necessarily to leading order $h \sim K/r$ for some constant $K$. Substituting this and the form for $f$  in the first equation I find $K^2 = (1+2C\lambda_c)/(4\lambda_h)$. The second equation then gives $C = (1+2\lambda_cK^2)/(2\lambda_f)$. There is a real solution $(C,K)$ as long as $\lambda_c^2 < 2\lambda_h\lambda_f$. If no mistakes, this suggests the ansatz $h \propto 1/r$ at large $r$ is inaccurate for larger $\lambda_c$. This likely means $h$ must not admit a power expansion at large $r$.
Nonetheless, numerically solving the system by imposing that $h\propto 1/r$ at large $r$ (to help the boundary value problem solver converge), I find the solution for $\lambda_c=\lambda_f=\lambda_h=1$ looks like this
 A: Your analysis is basically correct. If you insist on the asymptotic forms $f-1 \sim -C/r^2$ and $h\sim K/r$ at $r\to\infty$ with $K,C\ne 0$, then real values of $C$ and $K$ exist only when
$$ \frac{\lambda_f + \lambda_c}{2\lambda_f \lambda_h - \lambda_c^2} > 0. $$
Assuming $\lambda_{f,h,c} > 0$, reduces this to your inequality $\lambda_c^2 < 2\lambda_f\lambda_h$. So indeed no solution with these requirements exists for instance for large positive values of $\lambda_c$. The values of the parameters $\lambda_c = \lambda_f = \lambda_h = 1$ satisfy the inequality ($1 < 2$), so the desired asymptotic expansion exists and that's consistent with your numerical solution.
One thing you neglected to mention or didn't notice is that allowing $K=0$ fixes $C=\frac{1}{2\lambda_f}$ and corresponds to setting $h=0$ and $f$ to the solution of the second equation with $\lambda_c=0$. Provided this $f$ exists, setting $h=0$ both satisfies the full coupled system and the boundary conditions at $r=0$. This trick works for any value of $\lambda_c$, though it may be unsatisfactory for you for other reasons.
I suspect that trying to find a solution numerically with parameters say $\lambda_c=10$, $\lambda_f = \lambda_h = 1$, which fails to satisfy the inequality, will either fail, or converge to the "trivial" solution with $h=0$.
