How do you integrate an expression over a variable in the limit of an integral? I am trying to follow the steps to solve the integro-differential equation that arises from a plasma sheath problem given in this paper. This is the step I can't follow:

$$\epsilon_o\frac{d}{d\varphi}\biggl(\frac{E^2}{2}\biggr) = \sqrt{\frac{m_e}{2e}}\frac{j_{eo}}{\sqrt{\varphi}} - \sqrt{\frac{m_i}{2e}}\frac{j_{eo}}{\lambda_I}\int_\varphi^{\varphi_w}\frac{\varphi'/\varphi_I - 1}{\sqrt{\varphi' - \varphi}}\frac{d\varphi'}{E'},\tag{4}$$
where $\lambda_I =1/\sigma_o n_a$ is the ionization mean free path for the electron energy $E_I$. The integration of Eq. (4) over $\varphi$ leads to an integral equation for the electric field $E(\varphi)$,
$$\begin{multline}\frac{\epsilon_o}{4j_{eo}}\sqrt{\frac{e}{2m_e}}(E_w^2 - E^2) = \bigl(\sqrt{\varphi_w} - \sqrt{\varphi}\bigr)\\- \frac{1}{\lambda_I}\sqrt{\frac{m_i}{m_e}}\int_\varphi^{\varphi_w}\biggl(\frac{\varphi'}{\varphi_I} - 1\biggr)\sqrt{\varphi(z') - \varphi(z)}\frac{d\varphi'}{E(\varphi')},\tag{5}\end{multline}$$

The paper claims that instead of solving the integro-differential equation numerically from the form above, both sides of the equation can be integrated with respect to $\varphi$. I am not sure if it is valid to integrate both sides of this equation with respect to this variable since it appears in the lower limit of the integral. 
Can someone explain how you would handle integrating both sides when $\varphi$ is in the the limit of the integral? Thanks!
Or if anyone has an argument for why a mistake might have been made in this step that would be helpful too. 
 A: Consider an expression of the form
$$\int_x^c f(x, y)\mathrm{d}y$$
where $c$ is a constant. You can think of this as a mapping from real numbers to real numbers: you pick any real number $x$, plug it in, and calculate the value of the integral. That's exactly what a single-variable function is. So you can label this function $I_c$, define it as
$$I_c(x) = \int_x^c f(x, y)\mathrm{d}y$$
and then it should make sense that you can integrate it like any other function:
$$\int_a^b I_c(x)\mathrm{d}x = \int_a^b\int_x^c f(x, y)\mathrm{d}y\,\mathrm{d}x$$
So that tells you that integration over a variable that appears in the limit of an inner integral is a perfectly reasonable thing to do, and conceptually, there's nothing complicated about it.
Actually coming up with a symbolic expression for the double integral is another matter, of course. In this case, the major step in going from equations (4) and (5) is integrating over $\varphi$, so let's do that: for the term on the left side,
$$\int_{\color{red}\varphi}^{\varphi_w}\epsilon_o\frac{\mathrm{d}}{\mathrm{d}\varphi}\biggl(\frac{E^2}{2}\biggr)\mathrm{d}\varphi = \epsilon_o\biggl(\frac{E_w^2}{2} - \frac{\color{red}{E}^2}{2}\biggr)\tag{A1}$$ 
where $E_w = E(\varphi_w)$ and $\color{red}{E} = E(\color{red}{\varphi})$ and for the first term on the right,
$$\int_{\color{red}\varphi}^{\varphi_w}\sqrt{\frac{m_e}{2e}}\frac{j_{eo}}{\sqrt{\varphi}}\mathrm{d}\varphi = \sqrt{\frac{m_e}{2e}}j_{eo}\bigl(2\sqrt{\varphi_w} - 2\sqrt{\color{red}{\varphi}}\bigr)\tag{A2}$$ 
The reason I'm using some red variables, by the way, is that when a variable of integration appears as one of the limits, you should consider it to be a separate variable inside and outside the integral. So think of $\color{red}{\varphi}$ and $\varphi$ as different variables. If you prefer, you could give the variable a different label instead of using a different color (so you could call it, say, $\varphi_b$, instead of $\color{red}{\varphi}$).
Anyway, that was easy. The tricky term is
$$\sqrt{\frac{m_i}{2e}}\frac{j_{eo}}{\lambda_I}\int_{\color{red}{\varphi}}^{\varphi_w}\int_\varphi^{\varphi_w}\frac{\varphi'/\varphi_I - 1}{\sqrt{\varphi' - \varphi}}\frac{\mathrm{d}\varphi'}{E'}\mathrm{d}\varphi\tag{A3}$$
In this one you can't do the integral over $\varphi'$ because you don't know $E' = E(\varphi')$ and you don't have enough information about it to express it as something you can integrate.
What they've done in the paper is exchange the order of integration. This is a procedure you can use on any multiple integral where the limit of the inner integral depends on an outer variable of integration. For example, in a double integral of the form
$$\int_a^b \int_x^b f(x, y)\mathrm{d}y\,\mathrm{d}x$$
the region of integration is
$$\begin{align}
a &\leq x \leq b &
x &\leq y \leq b
\end{align}$$
which is the blue shaded region in this picture:

But you can express the same region as
$$\begin{align}
a &\leq y \leq b &
a &\leq x \leq y
\end{align}$$
as shown by this picture:

which gives you this identity
$$\int_a^b \int_x^b f(x, y)\mathrm{d}y\,\mathrm{d}x = \int_a^b \int_a^y f(x,y)\mathrm{d}x\,\mathrm{d}y$$
Using this procedure on equation (A3) from above turns it into
$$\sqrt{\frac{m_i}{2e}}\frac{j_{eo}}{\lambda_I}\int_{\color{red}{\varphi}}^{\varphi_w}\int_{\color{red}{\varphi}}^{\varphi'}\frac{\varphi'/\varphi_I - 1}{\sqrt{\varphi' - \varphi}}\frac{1}{E'}\mathrm{d}\varphi\,\mathrm{d}\varphi'$$
which is a fairly simple function of $\varphi$, integrable as $\int 1/\sqrt{\varphi'-\varphi}\mathrm{d}\varphi = -2\sqrt{\varphi'-\varphi}$. The full result after performing the inner integral over $\varphi$ is
$$\sqrt{\frac{m_i}{2e}}\frac{j_{eo}}{\lambda_I}\int_{\color{red}{\varphi}}^{\varphi_w}\biggl(\frac{\varphi'}{\varphi_I} - 1\biggr)\Bigl(-2\underbrace{\sqrt{\varphi' - \varphi'}}_{0} + 2\sqrt{\varphi' - \color{red}{\varphi}}\Bigr)\frac{\mathrm{d}\varphi'}{E'}\tag{A4}$$
Putting together all the terms, (A1), (A2), and (A4), we get
$$\begin{multline}\epsilon_o\biggl(\frac{E_w^2}{2} - \frac{\color{red}{E}^2}{2}\biggr) = 2\sqrt{\frac{m_e}{2e}}j_{eo}\bigl(\sqrt{\varphi_w} - \sqrt{\color{red}{\varphi}}\bigr) \\
- 2\sqrt{\frac{m_i}{2e}}\frac{j_{eo}}{\lambda_I}\int_{\color{red}{\varphi}}^{\varphi_w}\biggl(\frac{\varphi'}{\varphi_I} - 1\biggr)\sqrt{\varphi' - \color{red}{\varphi}}\frac{\mathrm{d}\varphi'}{E'}\end{multline}$$
Then divide both sides by $2\sqrt{\frac{m_e}{2e}}j_{eo}$ and you get equation (5), except for one factor of 2 on the left side which I can't seem to account for (maybe it's lost somewhere in my calculations).
A: Rather than deriving eq. (5) from eq. (4) via integration, it is enough to derive eq. (4) from eq. (5) via differentiation [and check that eq. (5) is satisfied when the integration limits coincide].
Hence differentiate on both sides of eq. (5) with the operator
$$ -2j_{eo}\sqrt{\frac{2m_e}{e}}\frac{d}{d\varphi}.$$
On the rhs. of eq. (5), one should differentiate under the integral sign. In principle, there is also a contribution from the lower limit of the integral, which turns out to be zero because the integrand vanishes in the lower limit. [Note that $\varphi(z^{\prime})$ in the integral is a typo. It should read  $\varphi^{\prime}$.] In the end eq. (4) emerges.
