$\color{green}{\textbf{Version of 15.08.20.}}$
At first,
$$I=\int_0^1x(1-x)f(x)f'(x)\text{d}x \,\overset{IBP}{=\!=\!=}\,\frac12x(1-x)f^2(x)\bigg|_0^1+\frac12\int_0^1 (2x-1)f^2(x)\text{ d}x,$$
$$I=\frac12\int_0^1 (2x-1)f^2(x)\text{ d}x.\tag1$$
Let
$$x = \frac{y+1}2,\quad \text{ d}x = \frac12\text{ d}y,\quad y = 2x-1,\quad g(y) = f\left(\frac{y+1}2\right), \tag2$$
then
$$f(x) = g(2x-1) = g(y),\quad f'''(x) = 8g'''(2x-1) = g'''(y),\tag3$$
$$I = \frac14\int\limits_{-1}^{1} y g^2(y)\text{ d}y,\tag4$$
under the conditions
$$g(-1) = g'(-1) = g(1) = 0,\quad |g'''(y)| \le \frac18.\tag5$$
Decomposition to the even and the odd parts
$$g(y)=g^\,_+(y)+g^\,_-(y),\quad g^\,_\pm(y) = \frac12(g(y)\pm g(-y)),\quad
g^\,_\pm(-y) = \pm g^\,_\pm(y),\tag6$$
gives
$$I = \int\limits_{0}^{1} y g^\,_+(y)\,g^\,_-(y)\text{ d}y.\tag7$$
In accordance with the Shwartz inequality,
$$I^2 \le \int\limits_{0}^{1} \big(y g^\,_+(y)\big)^2\text{ d}y\cdot \int\limits_{0}^{1} g^2_-(y)\text{ d}y,\tag8$$
wherein $(8)$ becames the equality if
$$|g^\,_-(y)| = y\,g^\,_+(y).$$
Then from $(5)$ should
$$g^\,_+(y) = (1-y^2)h(y),\quad \big|g^\,_-(y) \big| =(y-y^3)h(y)\tag{9}$$
Therefore, the function
$$g(y)=(1+y)(1-y^2) h(y)$$
maximizes $|I|$ under the conditions $(5)$ near $y=\pm1.$
Taking in account the rest of the conditions, one can get $h(y) = \text{constant} =\frac1{48},$
$$g(y) = \frac1{48}(1+y)(1-y^2),\tag{10}$$
$$48^2I_{opt} = \int\limits_0^1 (y^3-y)^2\text{ d}y
= \int\limits_0^1 (y^6-2y^4+y^2)\text{ d}y = \frac17-\frac25+\frac13 = \frac8{105},\tag{11}$$
$$\color{brown}{\mathbf{|I| \le \frac1{30240}}},$$
$$f_{opt}(x) = \pm g_{opt}(2x-1) = \pm \frac16 (x^2-x^3).$$