Calculus of variations with two functions and inequality

The proposed problem is typical for problems of the calculus of variations. A classical (exact) approach to solving similar problems requires us to present the integrand as a function of an independent variable, an optimized function, and its derivative. Abandonment of the use of the derivative allows us to consider the author's solution as a kind of trial and error method.

To use the classical approach, one can use the function $$G(y) = \int\limits_0^yg(y^*)dy^*,\qquad(1)$$ then $$\begin{cases} g(y) = G'(y)\\ Q(u) = \int\limits_0^h q(u)\,dy,\quad q(u) = u\\ B(u,G') = \int\limits_0^h b(u,G')\,dy,\quad b(u,G') = uG'\\ M(u,G) = \int\limits_0^h m(u,G)\,dy,\quad m(u,G) = u^2+G. \end{cases}\qquad(2)$$ Let us try to combine the elements of Lagrange multipliers method with theory of the Euler-Lagrange calculus of variations.

Lagrange function is $$F(u,G,\lambda,\mu) = Q + \lambda(B-B_0) + \mu(M-M_0),$$ the variations are $$\begin{cases} \dfrac{\delta q}{\delta u} = 1\\ \dfrac{\delta b}{\delta u} = G' = g\\ \dfrac{\delta b}{\delta G} = -\frac d{dy}u = -u'\\ \dfrac{\delta m}{\delta u} = 2u\\ \dfrac{\delta m}{\delta G} = 1, \end{cases}\qquad(3)$$ where have used Euler-Lagrange theorem of variations $$\dfrac{\delta F(y,f,f')}{\delta f} = F'_f - \dfrac{d}{dy}F'_{f'}$$ Equaling to zero $\dfrac{\delta F}{\delta u},\ \dfrac{\delta F}{\delta G},\ F'_\lambda,\ F'_\mu,\ $ one can get the system $$\begin{cases} 1+\lambda g + 2\mu u = 0\\ -\lambda u' + \mu = 0\\ \int\limits_0^hug\,dy = B_0\\ \int\limits_0^h(u^2 +G)\,dy = M_0. \end{cases}\qquad(4)$$ The second equation means $$\boxed{u = \frac{\mu}{\lambda}\,y+u_0},$$ then $$Q = \int\limits_0^h \left(\frac{\mu}{\lambda}\,y+u_0\right)\,dy = \frac{\mu }{2\lambda}h^2+u_0h,\qquad(5)$$

Parameters $B_0$ and $M_0$ can be used for calculation of unknowns $\lambda $ and $\mu$ from the system $(4).$

$$\begin{cases} g(y) = -2\frac\mu\lambda \left(\frac{\mu}{\lambda}\,y+u_0\right) - \frac1{\lambda},\\ G(y) = u_0^2-\left(\frac{\mu}{\lambda}\,y+u_0\right)^2 - \frac1{\lambda}y,\\ B_0 = \frac23\left(u_0^3 - \left(\frac{\mu}{\lambda}\,h+u_0\right)^3\right) + \frac1{2\mu}\left(u_0^2 - \left(\frac{\mu}{\lambda}\,h+u_0\right)^2\right),\\ M_0 = u_0^2h-\frac1{2\lambda}h^2, \end{cases}$$ $$\begin{cases} 4\left(\frac{\mu h}\lambda\right)^3 + 12u_0\left(\frac{\mu h}\lambda\right)^2 + 3\left(4u_0^2+\frac h\lambda\right)\left(\frac{\mu h}\lambda\right) +6u_0\frac h\lambda + 6B_0 = 0\\ \frac h\lambda = 2u_0^2 - 2\frac {M_0}h, \end{cases}$$ $$\left(\dfrac{\mu h}\lambda\right)^3 + 3u_0\left(\dfrac{\mu h}\lambda\right)^2 + \dfrac32\left(3u_0^2 - \dfrac {M_0}h\right)\left(\dfrac{\mu h}\lambda\right) + 3u_0\left(u_0^2 - \dfrac{M_0}h\right) + \dfrac32B_0 = 0$$ $$\left(\dfrac{\mu h}\lambda+u_0\right)^3 + \dfrac32\left(u_0^2 - \dfrac{M_0}h\right)\left(\dfrac{\mu h}\lambda + u_0\right) + \dfrac12\left(3B_0 - 3\dfrac {M_0}h u_0+u_0^3\right) = 0.\qquad(6)$$ The equation $(6)$ allows to obtain explicit expression for $\dfrac{\mu h}\lambda$ via $u_0$ - eg, using Vieta's substitution.

It's substitution in $(5)$ gives the requred function $Q(u_0)$.

UPD

Let us change the "fresh" conditions: $$ \begin{cases} \lambda>0\\ \mu\lambda<0\\ g(0)<g^* \end{cases}\rightarrow \begin{cases} u_0^2 > \dfrac{M_0}h,\\ \dfrac{\mu h}\lambda < 0\\ u_0 < \dfrac{\lambda g^* +1}{-2\mu}. \end{cases} $$ The first of them means that equation $(6)$ in the form of $$v^3+pv+q=0,$$ where $$v = \dfrac{\mu h}\lambda + u_0,\quad p = \frac32\left(u_0^2-\frac {M_0}h\right) > 0,\quad q = \dfrac12\left(3B_0 - 3\dfrac {M_0}h u_0+u_0^3\right),$$ have the single real root $$v = \sqrt[3]{-\dfrac q2+\sqrt{\dfrac{q^2}4 + \dfrac{p^3}{27}}} + \sqrt[3]{-\dfrac q2 - \sqrt{\dfrac{q^2}4 + \dfrac{p^3}{27}}},$$ аnd the reversal point of the sign of $\left(\dfrac{\mu h}\lambda\right)$ can be determined from the equation $$u_0^3 - \dfrac{M_0}h u_0 + \dfrac{B_0}2 = 0,$$ $$\dfrac23u_0p + B_0 = 0,$$ those the second condition should be performed for $u_0 > -\dfrac{3B_0}{2p}.$