I'm currently working on a nonlinear system of ODE's of the following form: $$f''(x) = F\left(f(x),f'(x),g(x),g'(x)g''(x)\right) $$ $$g''(x) = F(g(x),g'(x),f(x),f'(x),f''(x)) $$ with given initial boundary conditions for $f,g,f'$ and $g'$ at $x=0$. Since the system is quite complex and nonlinear, I'm pretty sure that a numerical method is needed in order to evaluate $f(x)$ and $g(x)$. For instance, the first equation is written as: $$f''\left\{ \left[A+\frac{Bg'}{\sqrt{1+g'^2}}\right] \left(\frac{g+2gf'^2}{\left(1+f'^2\right)^{3/2}}\right) - \frac{3Bff'\sqrt{1+g'^2}}{\left(1+f'^2\right)^{5/2}} \right\} = \left[A+\frac{Bf'}{\sqrt{1+f'^2}}\right]\sqrt{1+g'^2} + Cg - \\ B\left[\frac{f'\sqrt{1+g'^2}}{\left(1+f'^2\right)^{3/2}} + \frac{fg'g''}{\sqrt{1+g'^2}\left(1+f'^2\right)^{3/2}}+\frac{f'gg''}{\sqrt{1+f'^2}\left(1+g'^2\right)^{3/2}}\right]- \left[A+\frac{Bg'}{\sqrt{1+g'^2}}\right]\frac{gf'}{\sqrt{1+f'^2}} $$ where $A,B$ and $C$ are constants and the dependence on $x$ is omitted by simplicity. The second equation of the system is obtained simply by substituting all $f$ (and their derivatives) by $g$ (and their derivatives).
My doubt is on which numerical method should I adopt to solve the system? I have some experience with the RK4 method to solve second order ODEs, which is similar to a system of ODEs, but in this case, it looks different given that $f''(x)$ is a function of $g''(x)$ and vice-versa, which looks recursive to program. Any help provided will be very appreciated. Thank you in advance!
EDIT: The equations arise from the Euler-Lagrange equations of the following functional: $$J=\int_0^H{\left\{\left[A+\frac{Bf'}{\sqrt{1+f'^2}}\right]f\sqrt{1+g'^2} + \left[A+\frac{Bg'}{\sqrt{1+g'^2}}\right]g\sqrt{1+f'^2} + Cfg \right\}dx} $$