Solving a system of separable ODE's I am attempting to solve the system of first-order ODE's:
$$ \frac{dy_1(x)}{dx} = g(y_1)\frac{y_2(x)}{\int_{a}^b y_2(z) dz} $$
$$ \frac{dy_2(x)}{dx} = \frac{y_2(x)}{f(y_1(x),x)} \frac{\partial f(y_1(x),x)}{\partial x} $$
with boundary conditions $y_1(a) = c$, $y_1(b) = d$ for constants $a,b,c,d$.
My problem is conceptual, and very likely silly: I'm not sure whether I can break apart this problem apart in order to solve it. I would like to exploit separability to obtain a general solution for the second equation:
$$ y_2(x) = e^N f(y_1(x),x) $$
with N a constant, and use this to reduce the first equation to
$$ \frac{dy_1(x)}{dx} = g(y_1)\frac{f(y_1(x),x)}{\int_{a}^b f(y_1(z),z) dz} $$
which I would then solve numerically, presumably by treating the integral as a constant to be determined. I'm not seeing the problem with this approach, but on the other hand it strikes me as "too good to be true." My questions:

*

*Am I mistaken in thinking that $ y_2(x) = e^N f(y_1(x),x) $ is in fact the general solution for $y_2$, irrespective of the dependence of $y_1$ on $y_2$?

*When numerically solving a differential equation (say IVP) of the form $y_1'(x) = h\left(y_1,x,\int_{a}^b f(y_1(z),z) dz\right)$ with $a,b$ constant, is it problematic to treat the integral as a constant $A$, use an initial guess for $A$ + e.g Euler's method to solve for $y_1$, and then iterating until the guess and actual value of $A$ converge?

Any input would be much appreciated.
 A: Your system does not define $y_2(x)$ uniquely. If $(y_1(x), y_2(x))$ is a solution then also $(y_1(x), Cy_2(x))$ where $C \neq 0$ is a constant also is a solution.
Thus you may impose an extra condition on $y_2(x)$, the simplest to me is
$$
\int_a^b y_2(x) dx = 1.
$$
After that the system reduces to
$$
y_1'(x) = g(y_1(x)) y_2(x)\\
y_2'(x) = y_2(x) F(x, y_1(x))
$$
with conditions $y_1(a) = c, y_1(b) = d, \int_a^b y_2(x) dx = 1$.
Here I've denoted $F(x, y) = \frac{f_x(y, x)}{f(y, x)}$.
Having a system of two first-order ODE's with three conditions makes me think that the problem does not have a solution in general (but may have for some particular values of $c$ or $d$, like an eigenproblem).
The general solution to the second equation is
$$
y_2(x) = C_1 \exp\left(\int_a^x F(s, y_1(s)) ds\right)
$$
with $C_1$ determined by $\int_a^b y_2(x) dx = 1$ condition. Practically it is easier to put $y_2(a) = C_1 = 1$ and normalize $y_2(x)$ after integration.
Summarizing, the following algorithm may solve the problem:

*

*Take some initial $y_1(x)$, e.g. linear approximation $y_1(x) = c + \frac{x-a}{b-a} (d - c)$.

*Having $y_1(x)$ compute the ${\tilde y}_2(x)$ from the following Cauchy problem
$$
\begin{cases}
\tilde y_2'(x) = \tilde y_2(x) F(x, y_1(x)), &\quad x \in [a, b]\\
I'(x) = \tilde y_2(x), &\quad x \in [a, b]\\
\tilde y_2(a) = 1\\
I(a) = 0
\end{cases}
$$
Here $I(x) = \int_a^x y_2(s) ds$, so $I(b) = \int_a^b y_2(x) dx$.

*Compute $y_2(x)$ by $$y_2(x) = \frac{\tilde y_2(x)}{I(b)}.$$

*Plug obtained $y_2(x)$ into the first equation and solve the new Cauchy problem
$$
\begin{cases}
\tilde y_1'(x) = g(\tilde y_1(x)) y_2(x)\\
\tilde y_1(a) = c\\
\end{cases}
$$

*Now, for sure, $\tilde y_1(b) \neq d$, so we need to correct $y_1(x)$ and start from the item 2. And this is quite tricky. As I said, the problem almost certainly won't have a solution. If it was a boundary problem with three conditions and the functions we might have use the residual $\delta = d - \tilde y_1(b)$ to correct some free parameter (see shooting method), but there is no such parameter in this problem, at least as it was stated. You may try combining new approximation $y_1(x)$ as
$$y_1(x) := (1-\omega) y_1(x) + \omega \left[\tilde y_1(x) + \frac{x-a}{b-a}\delta\right]$$
and try to deduce how $\delta$ behaves in iterations depending on different $\omega$ values. The extra term in brackets is needed to correct the right boundary conditions in new $y_1(x)$ before restarting from item 2.

