Nonlinear Second-order ODE BVP with 4 boundary conditions

Question

My Lagrangian comes out in this form when I impose spherical symmetry:

$$ φ''(ρ)+{3\overρ} φ'(ρ)+{4μ^4\over M^2} φ(ρ)-{4μ^4\over M^4} φ^{3}(ρ)-{μ^4\over2M} ϵ=0 $$

The following boundary conditions apply: $ φ'(0)=0 $,$ φ'(∞)=0$, $ φ(0)=-φ_0$ And I've found that the fourth boundary value is:

$ φ(∞)={M(2^{1\over3}(\sqrt(81 \epsilon^2-768)-9 \epsilon)^{2\over3}+8 (3)^{1\over3})\over2 (6)^{2\over3} (\sqrt(81 \epsilon^2-768)-9 \epsilon)^{1\over3}}=φ_1$

The boundary on $ \phi(0) $ is a parameter I will find through variation according to specifics of my research, so for now I have simply left it as $\phi_0$. Assume everything is positive except the function, which can have negative values (thus the - sign in $φ(0)$). Also, $\rho$ is the only variable, so ODE.

I want to model the scalar field, but am having trouble solving for it. I know there is no analytic solution and I can easily redefine $φ'(ρ)=g(ρ)$ to get $2$ equations, but I am stumped about what numerical method to use to solve this monster.

How should I go about solving this?
Edit
Also, I wouldn't say no to tips on making this question more "answer-friendly".

Edit 2
I've found the following equation fits the boundary conditions, so that might be something to start with even if it doesn't fit the ODE:

$$\Phi(\rho)=\phi_1 e^{-\rho}+(2\phi_1+\phi_0)tanh({\phi_1 \rho\over 2\phi_1+\phi_0})-\phi_1-\phi_0$$

Answer 1

So, I've discovered I can use the Runge-Kutta 4th order method.

I define:

$$\phi'_{n+1}=\phi'_n+{1\over6}(K_1+2K_2+2K_3+K_4)$$ $$\phi_{n+1}=\phi_n+{1\over6}(L_1+2L_2+2L_3+L_4)$$ $$\rho_{n+1}=\rho_n+h$$ $K_1=({4\mu^4\over M^4}\phi_n^3-{4\mu^4\over M^2}\phi_n-{3\over\rho_n}\phi'_n+{\epsilon\mu^4\over2M})h$
$L_1=\phi'_nh$
$K_2=({4\mu^4\over M^4}(\phi_n+{1\over2}L_1)^3-{4\mu^4\over M^2}(\phi_n+{1\over2}L_1)-{3\over{\rho_n+{1\over2}h}}(\phi'_n+{1\over2}K_1)+{\epsilon\mu^4\over2M})h$
$L_2=(\phi'_n+{1\over2}K_1)h$
$K_3=({4\mu^4\over M^4}(\phi_n+{1\over2}L_2)^3-{4\mu^4\over M^2}(\phi_n+{1\over2}L_2)-{3\over{\rho_n+{1\over2}h}}(\phi'_n+{1\over2}K_2)+{\epsilon\mu^4\over2M})h$
$L_3=(\phi'_n+{1\over2}K_2)h$
$K_4=({4\mu^4\over M^4}(\phi_n+L_3)^3-{4\mu^4\over M^2}(\phi_n+L_3)-{3\over{\rho_n+h}}(\phi'_n+K_3)+{\epsilon\mu^4\over2M})h$
$L_4=(\phi'_n+K_3)h$

I have $\phi'_0=0$. Turns out, what I needed to do was vary my estimate for $\phi_0$ and keep running the simulation until I found an initial value that resulted in final values (as $\rho\rightarrow\infty$) that agree with the other two boundary conditions.