How do I solve this nonlinear ODE given the asymptotic series solutions as follows? Differential Equation:
$$-{\frac { \left( {\frac {\rm d}{{\rm d}R}}f \left( R \right)  \right) 
^{2}}{2\,f \left( R \right) }}+{\frac {{\rm d}^{2}}{{\rm d}{R}^{2}}}f
 \left( R \right) +{\frac {{\frac {\rm d}{{\rm d}R}}f \left( R
 \right) }{R}}+2\,f \left( R \right) -2\, \left( f \left( R \right) 
 \right) ^{2}-2\,{\frac {{B}^{2}}{f \left( R \right) {R}^{2}}}=0. $$
Series Solution at $R \to 0$:
$${R}^{-2}+1+ \left( {\frac {{B}^{2}}{3}}+{\frac{1}{3}} \right) {R}^{2}
+ \left( {\frac{2}{33}}+{\frac {2\,{B}^{2}}{33}} \right) {R}^{4}+
 \left( {\frac{31}{2277}}+{\frac {53\,{B}^{2}}{2277}}+{\frac {2\,{B}^{
4}}{207}} \right) {R}^{6}+ \left( {\frac{70}{29601}}+{\frac {136\,{B}^
{2}}{29601}}+{\frac {2\,{B}^{4}}{897}} \right) {R}^{8}+O \left( {R}^{
10} \right).$$
Series Solution as $R \to \infty$:
$$1-{\frac {{B}^{2}}{{R}^{2}}}+{\frac {-2\,{B}^{4}-2\,{B}^{2}}{{R}^{4}}}
+{\frac {-7\,{B}^{6}-23\,{B}^{4}-16\,{B}^{2}}{{R}^{6}}}+{\frac {-30\,{
B}^{8}-216\,{B}^{6}-474\,{B}^{4}-288\,{B}^{2}}{{R}^{8}}}+O \left( {R}^
{-10} \right) .$$
B is a free parameter. I tried solving it like a boundary value problem using series solution at very small R and large R as boundary conditions. I tried Newton iteration and imaginary time propagation for that. But that worked only for B=1.
If I try solving it as an initial value problem starting from some Rmax, even RKF45 doesn't give good result.
Can someone please suggest either analytical or numerical way of solving this equation? Or perhaps a way of analyzing the properties of this differential equation other than frobenius series solution?
(If the questions lacks details, please let me know before downvoting)
Edit: 2 downvotes without any explanation. I mean if you don’t wanna respond, then don’t respond. Why do you have to ruin my chances of getting any help? This is the worst forum. Most of the times people just keep downvoting without any explanation.

First, I'd prefer do a simple transformation $f(R) = (u(R))^2$ suggested by @EliBartlett.
After some manipulation with computer algebra system I got
$$
R\frac{d}{dR}\left(R \frac{du}{dR}\right) = R^2 (u^3 - u) + \frac{B^2}{u^3}
$$
Frobenius analysis shows that $u \sim \frac{1}{R}$ near $R = 0$. Let's introduce $x = \log R$ and $w = \log u$. The domain for $x$ becomes $(-\infty, -\infty)$ and the boundary conditions $w \to -x$ when $x \to -\infty$ and $w(+\infty) = 0$. Also $R \frac{d}{dR} = \frac{d}{dx}$.
The equation becomes
$$
\frac{d^2}{dx^2} e^w = e^{2x} (e^{3w} - e^w) + B^2 e^{-3w}.
$$
Consider a finite interval $x \in [-L, L]$ and a regular grid with step $h = \frac{2L}{N}$. The discrete problem becomes
$$
\frac{e^{w_{n+1}} - 2e^{w_n} + e^{w_{n-1}}}{h^2} = 
e^{2x_n} (e^{3w_n} - e^{w_n}) + B^2 e^{-3w_n}, \quad n = 1, \dots, N-1.\\
w_0 = -x_0, \quad w_N = 0.
$$
To avoid numerical cancellation let's divide the $n$-th equation by $e^{w_n}$:
$$
\frac{e^{w_{n+1} - w_n} - 2 + e^{w_{n-1} - w_n}}{h^2} = 
e^{2x_n} (e^{2w_n} - 1) + B^2 e^{-4w_n}, \quad n = 1, \dots, N-1.\\
w_0 = -x_0, \quad w_N = 0.
$$
This is a system that can be plugged into the Newton solver.
I've implemented this in Python (can be found here). Unfortunately, it works only when $B \lesssim 1.117$. I tried to make Newton's method robust by introducing safety factors, but no luck. I also tried scipy nonlinear solvers and they also hit a wall of no-convergence near $B = 1.117$. I feel that the problem has no solution for such $B$.
