Find $f(x)$ such that it maximizes $\int_0^1 \left(f(f(x))-f(x)\right) dx$ Inspired by this post Find the maximun of the sum $\sum_{k=1}^{n}(f(f(k))-f(k))$,
I came to the below question.

Find the $f(x)$ that satisfies:
$1.$ $0\le f(x)\le 1$
$2.$ Increasing for $\left(0,1\right)$
$3.$ Maximizes  $\quad\int_0^1 \left(f(f(x))-f(x)\right) dx$

After I tried a few elementary functions, now I am feeling that there might not be an analytic solution for the $f(x)$ (but I'm not sure). If that's the case, can we at least know the maximum of the integral?

I posted this question because the original discrete version of the question didn't seem to have a clear solution I was wondering if this continuous version of the question has any analytic solution. The accepted answer suggests that the solution could be the piecewise linear function. However, there's no clear proof or reasoning why that should be the optimal case.

What I've tried was:
$1$. I started from $f(x)=\frac{x+1}2$ as it was assumed in the original question I likned above.
$2$. I tried a few elementary functions such as $f(x)=x^p$, and then $f(x)=m\ln\left(\frac{x+a}{1+a}\right)+1$, and two terms, etc. (assumed $f(1)$ should be $1$ to maximize the integral) and could see that the integral kept increasing for more variables, so I thought there must be an optimum (but differentiable) $f(x)$ and that it might not be an analytically obtained due to $f(f(x))$. At that time I didn't know that I was only approaching to the step function.

Thanks for all the comments and the answers. Now I can see that the solution is probably the piecewise step function in the answer, so this question may not be as interesting as I thought. However, I think that we are still missing the proof that it is the optimal $f(x)$.
 A: This isn't an answer but I have a function which gets you as close to $\frac{1}{4}$ as you desire. Further, I think the method I used could be used to find an even higher value (EDIT: maybe not, see comments at bottom) but I'm not sure (I might come back to this if I can get better results). To this end, I will first give my example and then explain how I found it in the hopes that it might lead to a final answer at some point.
Define the function
$$f(x)=\begin{cases}
 \epsilon +\frac{1}{2} & x\leq \frac{1}{2}-\epsilon  \\
 -\frac{x^3 (1-2 \epsilon )}{\epsilon ^3}+\frac{3 x^2 (\epsilon -1) (2 \epsilon -1)}{2 \epsilon ^3}-\frac{3 x \left(4 \epsilon ^2-4 \epsilon +1\right)}{4 \epsilon ^3}-\frac{-8 \epsilon ^3-6 \epsilon ^2+5
   \epsilon -1}{8 \epsilon ^3} & \frac{1}{2}>x>\frac{1}{2}-\epsilon  \\
 1 & x\geq \frac{1}{2}
\end{cases}$$
for some $0<\epsilon<\frac{1}{2}$. Now, it is rather difficult to show but this function fits the conditions presented above. Further, we have
$$\int_0^1(f(f(x))-f(x))dx=\frac{1}{4}-\frac{3}{4}\epsilon+\frac{1}{2}\epsilon^2 $$
Clearly, we can get within any $\epsilon>0$ of $\frac{1}{4}$ as we desire.
Now, my method to find this function:
$1)$ I started under the assumption that $f(x)$ was a piecewise function
$$f(x)=\begin{cases}
 a_1 & x< b \\
 a_2 & x>b
\end{cases}$$
(ignoring what happens at $x=b$).
$2)$ I then plugged this function into mathematica and got the result
$$\int_0^1(f(f(x))-f(x))dx=(a_2-a_1)b$$
when $b<a_1$. For other possible cases, I got a negative value and since $f(x)=0$ gives an integral of $0$ I knew that these cases could be discarded.
$3)$ I looked at this and thought about how to maximize it. Obviously, $a_2$ should be set to $1$. This then gave me
$$=(1-a_1)b$$
when $a_1>b$.
$4)$ I then set $a_1=b+\epsilon$ to get
$$=(1-b-\epsilon)b$$
This function is maximized at $b\approx \frac{1}{2}$. Thus, our original function becomes
$$f(x)=\begin{cases}
 \frac{1}{2}+\epsilon & x< \frac{1}{2} \\
 1 & x>\frac{1}{2}
\end{cases}$$
$5)$ To make our function fit, I used the cubic
$$h(x)=sx^3+tx^2+ux+v$$
and set
$$f(x)=\begin{cases}
 \epsilon +\frac{1}{2} & x\leq \frac{1}{2}-\epsilon  \\
 h(x) & \frac{1}{2}>x>\frac{1}{2}-\epsilon  \\
 1 & x\geq \frac{1}{2}
\end{cases}$$
All that was left was to solve for $s,t,u,v$ such that $f(x)$ met all the conditions above. I then plugged $f(x)$ back into mathematica to get the answer above.
Commentary: The obvious next step is to try this process with a step function of three levels rather than two. However, this already seems extremely difficult and I fear some clever idea is needed in order to make this work. The hope eventually would be to extend this to a step function of $n$ steps, and then in the limit we would hopefully have a function which maximizes the integral.

EDIT: I tried to extend the method above to a step function with three steps, but I kept running up on the $\frac{1}{4}$ barrier. Don't know if that means it is impossible to get past that, but I'm beginning to think that may be the upper limit. Especially since the linked post has answers which have $\frac{n^2}{4}$ in them as a sort of limiting factor. That is in the limit
$$\lim_{n\to\infty}\frac{n^2}{4}\cdot \frac{1}{g(n)}=1$$
and
$$g(n)<\frac{n^2}{4}$$
(where in that question $g(n)$ was the maximum). This is kind of like the answer that I got.

EDIT 2: After condition $2$ was changed to monotonicly increasing, there exists a function $f(x)$ such that the integral is exactly $\frac{1}{4}$. For
$$f(x)=\begin{cases}
 \frac{1}{2} & x< \frac{1}{2} \\
 1 & x\geq \frac{1}{2}
\end{cases}$$
we have
$$\int_0^1(f(f(x))-f(x))dx=\frac{1}{4}$$
