I have a problem of the following structure: For some scalar $g$, functions $F(z)$ and $h(z)$ defined on $[0,\bar{z}]$ , and a non-linear operator $\phi(F,z)$ (in reality, $F$ and $h$ are vector valued)
The $F(z)$ and $g$ must fulfill the ODE st. boundary values: $$ 0 = -F'(z;g) + \phi(F,z;g)\\ F(0) = 0\\ F(\bar{z}) = 1 $$ and the equilibrium condition: $$ 0 = -c + \int_{0}^{\bar{z}}h(z;g)F'(z;g) d z $$
My question is: what the appropriate algorithm to solve this numerically? I can solve the non-linear ODE with finite differences (and some shooting method in a more complicated variation), but a nested iteration with the integral equilibrium condition is difficult and expensive.
What about collocation? For example, pick the $N$th order chebyshev approximation. The chebyshev roots and polynomials (adapted to the $[0,\bar{z}]$ bounds) are $\{z_n\}_{n=0}^{N}$ such that $z_0 = 0$ and $z_N = \bar{z}$ and the polynomial of order $k$ evaluated at the $n$ root is $T_{nk} = T_n(z_n)$ and the vector of all orders at the $n$th root is $T_n$
Then if I assume a vector $a_0,\ldots a_N$ and a vector of these $F(z) \approx \sum_{a=0}^{N} a_n T_n(z)$ and therefore, $$F(z_n) = a \cdot T_n$$
Since chebyshevs are easy to differentiate, assume another set of values such that $F'(z_n) = a \cdot U_n$
Then I can apply the ODE and boundary values to get the usual system in $a$. For the integral constraint, I already have the chebyshev nodes, so if I assume weights $w_k$ for the usual chebyshev quadrature, then I get a final equation in $g$ $$ 0 = -c + \sum_{k=0}^N w_k h(z_k;g) \, a \cdot U_k $$
Does that sound like a defensible algorithm (using a nonlinear solver to find the $a,g$)? Are there other approaches which might work better?