I'm assuming you are interested in the region $R<1$; the region $R>1$ as shown in the graph.
Dominique is right, you do need a second constraint or boundary condition. More about that below.
Both the wicked behavior near 0 and the combination $f''+\frac{2}{R}f'$ argue for making the substitution $x = \frac{1}{R}$. This transforms the equation to
$$
\frac{d^2f}{dx^2} = \frac{K_+}{x^3}\left( \sqrt{f} - \lambda_+\sqrt{1-f} \right)
$$
where per the notation on the graph, $K_+$ is $0.7$ and $\lambda_+$ is $0.3$.
This equation is to be solved numerically for $x > 1$. In that region, the equation is quite well behaved, and Runge-Kutta will work well. Depending on how near to zero you need to get, you will have to integrate to some value $x_f > 1$. In order to get the solution down to $R=\epsilon$ you would integrate to $x_f = 1/\epsilon$.
And now the zinger: What is that second boundary condition?
If you are given $\frac{df}{dR}$ at $R=1$, this is an initial value problem starting from $x=1$, and that is straightforward. For example, the dashed curve on the graph looks to have $$\left. \frac{df}{dR} \right|_{R=1} = K_+ = 0.7$$
However, the main interesting feature on the graph is that $f(R=0.3) = 0$. Your graph contains multiple distinct dashed curves (which might be one reason you found it so hard to reproduce) since the dashed curve to the left of $0.3$ does not match the one to the right. This suggests that the boundary conditions used were
$$
\left\{
\begin{array}{c}
f(1) = 1 \\ f(\lambda_+) = 0
\end{array}
\right.$$
Runge-Kutta is not ideal for a BVP but you can do an iterative solution by shooting for the fixed point at $x=1/0.3, f(x) = 0$, or there are relaxation methods that will work fine because the function is well behaved in this region.