Numerically Solving a Second Order Nonlinear ODE Okay, I have this not so pretty 2nd order non-linear ODE I should be able to solve numerically. 
$$f''(R) + \frac{2}{R} f'(R)=\frac{0.7}{R} \left( \frac{1}{\sqrt{f(R)}} - \frac{0.3}{\sqrt{1-f(R)}} \right),$$
$$f(1)=1.$$
The function around the origin is behaving very wildly.
I was thinking of breaking this guy up into a system of two first order ODE's and then solve, but I have no idea how to set this up. What method should I use to set up the system of ODE's?
If there is some other method rather than numerically solving a system of differential equations, please feel welcome to share. Thanks.

 A: The general method to reinterpret a higher-order ODE as a system of first order ODEs is to regard the derivatives of the unknown function as additional unknown functions. In your case, regard $f'$ as a new function $g$. Then the system of first order ODEs consists of two equations, the first being the original equation with $f'$ replaced by $g$ and $f''$ replaced by $g'$, and the second equation being $f'=g$.
A: 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. 
