I am looking for ways to solve the following system of boundary value implicit ODEs over the interval $t \in [0, 1]$: \begin{equation} \lambda (Fg - Gf)^3 + 4 FGfg(g-f) = 0 \\ fg(Fg-Gf) + 2FG(gf' - fg') = 0, \end{equation} where $f = F'$ and $g=G'$ and $f, -g \geq 0$. The boundary values are $F(0) = G(1) = 0$ and $F(1) = G(0) = 1$.

The equations arise as the Euler Lagrange equations of the following constrained optimization problem: \begin{equation} \text{minimize}_{F, G} \int_0^1 \frac{dx}{F(x) / f(x) - G(x) / g(x)}, \end{equation} subject to: \begin{equation}\frac{1}{2} \int_0^1 [f(x) \log \frac{2f(x)}{f(x) - g(x)} - g(x) \log \frac{2g(x)}{g(x) - f(x)}] dx = J_0 \end{equation}

The constraint is simply the Jensen Shannon divergence of the two distributions on $[0, 1]$, whose densities are given by $f$ and $-g$ respectively. The objective function is another similarity function between the two distributions related to the Jaccard distance. Note that only the minimization problem above is nontrivial. The maximization problem can be solved with a pair of distributions on a 2-point space. Finite dimensional numerical solutions suggest the minimizing $F$ and $G$ are strictly monotonic over $[0, 1]$, thus describe distributions on an infinite space.

Any analytic insight (not necessarily leading to a solution) would be greatly appreciated. Alternatively, suggestion on what commercial or open source solvers to use would be just as awesome: I have looked into implicit Runge-Kutta implementation in R and matlab but they seem pretty limited in scope, in particular, only capable of solving index-1 initial value problems, unless there is a way to convert BVP to IVP that I missed?