I wonder if the following equation could be further simplified by Mathematica. There are $22$ Integrals involved and many of them are in the range $y_l$ and $y_u$. I would guess that at least those terms could be collected in the same integral. I used the command FullSimplify but the result is the same with the input. Are there some more clever ways to simplify such type of equations?
(2*(1 - \[Epsilon]0)^2*Integrate[((Sqrt[f1[y]*l[yu]] - Sqrt[f0[y]*l[yl]
*l[yu]])*(-Sqrt[f1[y]*l[yl]] + Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] +
Sqrt[l[yu]])^2, {y, yl, yu}] - 2*(Integrate[(Sqrt[f0[y]]*(-Sqrt[f1[y]*l[yl]] +
Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] + Sqrt[l[yu]]), {y, yl, yu}] +
Integrate[f0[y], {y, -Infinity, yl}]*Sqrt[l[yl]])*(Integrate[(Sqrt[f0[y]]*
(Sqrt[f1[y]*l[yu]] - Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] + Sqrt[l[yu]]),
{y, yl, yu}] + Integrate[f0[y], {y, yu, Infinity}]*Sqrt[l[yu]]) -
Sqrt[(-2*(1 - \[Epsilon]0)^2*Integrate[((Sqrt[f1[y]*l[yu]] -
Sqrt[f0[y]*l[yl]*l[yu]])*(-Sqrt[f1[y]*l[yl]] +
Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] + Sqrt[l[yu]])^2, {y, yl, yu}] +
2*(Integrate[(Sqrt[f0[y]]*(-Sqrt[f1[y]*l[yl]] +
Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] + Sqrt[l[yu]]), {y, yl, yu}] +
Integrate[f0[y], {y, -Infinity, yl}]*Sqrt[l[yl]])*(Integrate[(Sqrt[f0[y]]*
(Sqrt[f1[y]*l[yu]] - Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] +
Sqrt[l[yu]]), {y, yl, yu}] + Integrate[f0[y], {y, yu, Infinity}]*
Sqrt[l[yu]]))^2 - 4*((Integrate[(Sqrt[f0[y]]*(-Sqrt[f1[y]*l[yl]] +
Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] + Sqrt[l[yu]]), {y, yl, yu}] +
Integrate[f0[y], {y, -Infinity, yl}]*Sqrt[l[yl]])^2 -
(1 -\[Epsilon]0)^2*(Integrate[(Sqrt[f0[y]]*(-Sqrt[f1[y]*l[yl]] +
Sqrt[f0[y]*l[yl]*l[yu]])^2)/(-Sqrt[l[yl]] + Sqrt[l[yu]])^2, {y, yl, yu}] +
Integrate[f0[y], {y, -Infinity, yl}]*l[yl]))*((Integrate[(Sqrt[f0[y]]*
(Sqrt[f1[y]*l[yu]] - Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] +
Sqrt[l[yu]]), {y, yl, yu}] + Integrate[f0[y], {y, yu, Infinity}]*
Sqrt[l[yu]])^2 - (1 - \[Epsilon]0)^2*(Integrate[(Sqrt[f1[y]*l[yu]] -
Sqrt[f0[y]*l[yl]*l[yu]])^2/(-Sqrt[l[yl]] + Sqrt[l[yu]])^2, {y, yl, yu}] +
Integrate[f0[y], {y, yu, Infinity}]*l[yu]))])/(2*((Integrate[(Sqrt[f0[y]]*
(Sqrt[f1[y]*l[yu]] - Sqrt[f0[y]*l[yl]*l[yu]]))/(-Sqrt[l[yl]] + Sqrt[l[yu]]),
{y, yl, yu}] + Integrate[f0[y], {y, yu, Infinity}]*Sqrt[l[yu]])^2 -
(1 - \[Epsilon]0)^2*(Integrate[(Sqrt[f1[y]*l[yu]] -
Sqrt[f0[y]*l[yl]*l[yu]])^2/(-Sqrt[l[yl]] + Sqrt[l[yu]])^2, {y, yl, yu}]
+Integrate[f0[y], {y, yu, Infinity}]*l[yu])))