Disclaimer: I'm a physicist, so I'm not sure if this kind of question is right for this board, or physics, or mathematica, or even overflow, and I don't know what exactly to search for, so hopefully someone can give me some guidance, or move the thread someplace else.
To get straight to the point: I have the following system of first order nonlinear ODEs:
$$h'=\frac{g}{1-\frac{R_h}{r}}$$ $$g'=-\frac{2}{r}g+\frac{dV[h(r)]}{dh(r)}$$ $$m'=4\pi r^2\bigg[\frac{1}{2}\frac{g^2}{1-\frac{R_h}{r}}+V[h(r)]\bigg]$$
Here $h,\, g,\, m$ are all functions of $r$, and $V$ is a known function of $h(r)$, while the boundary conditions are:
$$m(R_h+\epsilon)=0$$ $$g(R_h+\epsilon)=\frac{dV[h(r)]}{dh(r)}\bigg|_{R_h}\epsilon$$ $$h(R_h+\epsilon)=h(R_h)+\frac{R_h}{1-8\pi\alpha R_h^2 V[h(R_h)]}\frac{dV[h(r)]}{dh(r)}\bigg|_{R_h}\epsilon$$
Here $\epsilon\ll R_h$ (known), and $\alpha$ is a known constant. The problem is, the value of $h(R_h)$ is not known, but the (physical) requirement is that $\lim\limits_{r\rightarrow \infty}h(r)=0\quad (1)$. I've tried rewriting the equations in terms of $x=\frac{1}{r}$, but I get issues with infinities. So, what I did was solve them numerically for a bunch of initial values of $h(R_h)\equiv h_h$, and see what they look like. I get the following result: h(r) as a function of r for $R_h=0.1$
As can be seen, there seems to be a critical value of $h_h$ (approx. the yellow curve) for which condition $(1)$ is satisfied. What I'm looking for is a way to numerically find this value, and up to arbitrary accuracy as it's too time consuming to do it by hand for any value of $R_h$. I had a couple of ideas in mind, but I'm not entirely certain if I'm on the right track.
1) as the figure suggests, there seems to be a point where the solution goes from oscillating around +1 to oscillating around -1, so I was thinking of implementing a binary search algorithm, though I'm not sure which condition I should set to pick an interval (should I choose $\langle \frac{a+b}{2},b\rangle$ or $\langle a, \frac{a+b}{2}\rangle$)?
2) a different approach would be as follows: the figure suggests that the value that must be minimized is the area under the curve of $h(r)^2$. Hence, we can demand that $I=\int\limits_a^b h(r)^2\, dr$ is minimized, but that isn't the whole story as $h$ must satisfy the original differential equations, so I'm guessing Lagrange multipliers must be introduced into the mix somewhere, but I'm not sure how.
I'm currently more inclined to use option 1), but some help in the conditions would be appreciated. Of course, alternative ideas for finding the critical value of the boundary condition are welcome.