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?
Appropriate algorithm for (non-linear) ODE with integral equilibrium constraint: collocation?
algorithms;ode