Maximizing the value of $\int_0^1 f(x)f^{-1}(x)\ \mathrm dx$
Prove $\int_0^1 f(x)f^{-1}(x)dx \leq \frac{1}{3}$.
First we notice that $f(x) \leq x$ leads to $x \leq f^{-1}(x)$, and $f(x) \geq x$ leads to $x \geq f^{-1}(x)$, so $[f(x)-x][f^{-1}(x)-x] \leq 0$.
Integrate it.Then we get $\int_0^1 f(x)f^{-1}(x)dx+\int_0^1 x^2 dx \leq \int_0^1 x(f(x)+f^{-1}(x))dx=\int_0^1 xf(x) dx+\int_0^1 xf^{-1}(x)dx$.
Let $y=f^{-1}(x)$,then the second integral in the right is$\int_0^1 yf(y)df(y)=\frac{1}{2}\int_0^1 ydf^2(y)=\frac{1}{2}[1-\int_0^1 f^2(y)dy]$.
So$\int_0^1 f(x)f^{-1}(x)dx+\int_0^1 x^2 dx \leq \frac{1}{2}[1+\int_0^1 f(x)(2x-f(x))dx] \leq \frac{1}{2}[1+\int_0^1 x^2 dx]$.
done.
"=" iff $2x-f(x)=f(x)$,i.e. $f(x)=x$
The answer is indeed 1/3, which can be proved using the Fenchel-Young inequality for Legendre transforms.
Define $F(t):=\int_0^t f(x)dx$ so that $F$ is convex on $[0,1]$. The Legendre transform of $F$ is given by $G(t)=\sup_{u\in [0,1]} (tu-F(u))=\int_0^t f^{-1}(x)dx$ for $t \in [0,1]$.
Young's inequality (also called Fenchel's inequality) says that $ab \leq F(a)+G(b)$ for any $a,b \in [0,1]$.
Consequently we see that $f(x)f^{-1}(x)\leq F(f(x))+G(f^{-1}(x))$. Now notice from Fubini that $$\int_0^1F(f(x))dx = \int_0^1\int_0^{f(x)}f(u)dudx $$$$= \int_0^1 \int_0^1 f(u)1_{\{u<f(x)\}}dxdu=\int_0^1 f(u)(1-f^{-1}(u))du,$$ and symmetrically we obtain that $$\int_0^1G(f^{-1}(x))dx = \int_0^1 f^{-1}(u)(1-f(u))du.$$ Now integrating the identity $f(x)f^{-1}(x)\leq F(f(x))+G(f^{-1}(x))$ from $0$ to $1$, we get that $$\int_0^1 f(x)f^{-1}(x)dx \leq \int_0^1 f(u)(1-f^{-1}(u))du+\int_0^1 f^{-1}(u)(1-f(u))du.$$ Since $\int_0^1f(u)du+\int_0^1 f^{-1}(u)du=1$, the previous expression reduces to $$3\int_0^1 f(x)f^{-1}(x)dx\leq 1.$$
I believe you can prove it using the calculus of variations. Here is a sketch of how it would go.
Consider $h(x) = f(x) + \delta g(x)$ where $g(x)$ is a function in some class of nicely-behaved functions and $\delta$ is a small number. Then $h^{-1}(x) \approx f^{-1}(x) - \delta g(x)$. Now consider
$$u(\delta) := \int_0^1 h(x)h^{-1}(x) = \int_0^1 (f(x) + \delta g(x)) (f^{-1}(x) - \delta g(x))$$
Then
$$u'(\delta)|_{\delta = 0} = \int_0^1 (g(x)f^{-1}(x) - g(x) f(x)) = \int_0^1 g(x)(f(x) - f^{-1}(x))dx$$
Since this must be zero for all $g$, you get $f^{-1}(x) = f(x)$. But since $f$ is increasing, you get
$$x \le f(x) \le f(f(x)) = x$$
and therefore $f(x) = x$, so your answer is correct.