Se puede resolver con el cálculo de variaciones. Vamos a definir un funcional:
$$I(f) = \int_0^1 f(x)-f(x^2)^2 \ dx.$$
Entonces:
$$I(f+h) = \int_0^1 f(x)-f(x^2)^2 + h(x)-h(x^2)^2 -2f(x^2)h(x^2) \ dx = I(f) + \int_0^1 h(x)-2f(x^2)h(x^2) \ dx + o(h),$$
así que, con un cambio de variables,
$$\delta I (f) (h) = \int_0^1 h(x)-2f(x^2)h(x^2) \ dx = \int_0^1 h(x) \left(1-\frac{f(x)}{\sqrt{x}}\right) \ dx.$$
Obviamente, $f(x) = \sqrt{x}$ es un punto crítico, y resulta que $I(\sqrt{x}) = 1/3$. A ver que es un máximo, es suficiente con mirar la segunda derivada, o simplemente para factorizar toda la qudratic forma $I$:
$$I(\sqrt{x}+h(x)) = \frac{1}{3} -\int_0^1 h(x^2)^2 \ dx,$$
de manera continua $f(x) = \sqrt{x}+h(x)$, $I(f) = 1/3$ si y sólo si $h = 0$, por lo que si y sólo si $f(x)=\sqrt{x}$.