The answer is yes. First, to prove that $F$ has exactly one local maximum in $(0,\infty)$, it is enough to show that $F$ is increasing-decreasing on $(0,\infty)$ -- that is, for some $c\in[0,\infty]$ the function $F$ is non-decreasing on $(0,c]$ and non-increasing on $[c,\infty)$. $\big($Indeed, $F$ is clearly real-analytic on $(0,\infty)$ and hence can be constant on any non-empty open subinterval of $(0,\infty)$ only if it is constant on $(0,\infty)$. But $F(0+)=0$, whereas $F>0$ on $(0,\infty)$.$\big)$
Thus, it is enough to show that $F$ is the pointwise limit of a sequence of increasing-decreasing functions on $(0,\infty)$.
Let $g:=f^{-1}\colon[0,b]\to[0,1]$, where $b:=f(1)\in(0,\infty)$. Then $g$ is a convex increasing function such that $g(0)=0$. Let $g'$ be the left derivative of $g$, with $g'(0):=0$. Then $g'\colon[0,b]\to\mathbb{R}$ is a non-decreasing function, which is left-continuous on $(0,b]$. So, for all $t\in[0,b]$
$$g(t)=\int_0^t g'(u)\,du=\int_0^t\,du\int_{[0,u)}\,dg'(a)=\int_{[0,t)}\,dg'(a)\int_a^t du
=\int_{[0,b)}\,dg'(a)\,(t-a)_+,
$$
by Fubini's theorem, where $x_+:=0\vee x$.
Hence, by the change of variables $t=f(x)\iff x=g(t)$ and again Fubini's theorem,
$$
y^{-3/2}F(y) =\int_{[0,b]} \frac{dg(t)}{(y+t)^2}= \int_{[0,b)}dg'(a)\, \int_{[0,b]} \frac{d_t(t-a)_+}{(y+t)^2}
=\int_{[0,b)}dg'(a)\, \int_a^b \frac{dt}{(y+t)^2}
$$
$$
=\int_{[0,b)}dg'(a)\,\frac{b-a}{(y+b)(y+a)}
$$
for $y\in(0,\infty)$.
So, it is enough to show that any mixture $\big($that is, any finite linear combination with nonnegative coefficients$\big)$ of the functions of the form
$$(*)\qquad(0,\infty)\ni y\mapsto\dfrac{y^{3/2}}{(y+b)(y+a)}\quad\text{with a fixed $b\in(0,\infty)$ and all $a\in(0,b)$}
$$
is increasing-decreasing.
By rescaling $z:=y/b$, without loss of generality (wlog) $b=1$.
Now the crucial ingredient, Theorem 4 in Kemperman-modes, tells us that it suffices to consider just the mixtures of two functions of the form $(*)$, with $b=1$. Such a function is wlog given by the formula
$$ h(y)=\dfrac{y^{3/2}}{(y+1)(y+a_1)}+w\,\dfrac{y^{3/2}}{(y+1)(y+a_2)}$$
for some $a_1$ and $a_2$ in $(0,1)$, some $w>0$, and all $y>0$.
Note that
$$
h_1(y):=h'(y)\,\frac{2 (y+1)^2 \left(a_1+y\right){}^2 \left(a_2+y\right){}^2}{\sqrt{y}}
=a_1 \left(2 a_2 (w+1) (y+3) y+a_2^2 (y+3)+y^2 (-2 w y+2 w+y+3)\right)+y
\left(a_2 y (w (y+3)-2 y+2)-a_2^2 (y-1)-(w+1) (y-1) y^2\right)+a_1^2 w
\left(a_2 (y+3)-(y-1) y\right)
$$
equals $h'(y)$ in sign for all $y>0$. So, it is enough to show that there are no such $y_1$ and $y_2$ that $0<y_1<y_2$ and $h_1(y_1)<0<h_1(y_2)$. Since $h_1$ is algebraic (even a polynomial), this can be done completely algorithmically, say using a theory by Tarski. In Mathematica, this algorithm is implemented into Reduce[] and similar commands. As shown in the Mathematica notebook at Reduce.Mathematica (the corresponding pdf file at Reduce.pdf), the Reduce command reduces the conjunction of the conditions $0<y_1<y_2$, $h_1(y_1)<0<h_1(y_2)$, $0<a_1<1$, $0<a_2<1$, and $w>0$ to False, in about 18 sec.
This contradiction completes the proof.
Addendum: I got curious as to how essential the exponent $3/2$ in the definition of the function $F$ in question is. It turns out that the unimodality result holds if the exponent $3/2$ is replaced by any $p\in[1,2)$, and only in that case. Indeed, the functions
$$(**)\qquad(0,\infty)\ni y\mapsto g_{p,a}(y):=\dfrac{y^p}{(y+1)(y+a)}\quad\text{with $a\in[0,1]$}
$$
(cf. $(*)$)
are clearly monotonically decreasing if $p\in(-\infty,0]$ and monotonically increasing if $p\in[2,\infty)$. It remains to consider $p\in(0,2)$. Assume in the sequel that $(p,a,y)\in(0,2)\times[0,1]\times(0,\infty)$.
Then, for each pair $(p,a)$, there is a unique (mode) $m_{p,a}\in(0,\infty)$ such that
$$(0)\qquad g'_{p,a}(y)
\begin{cases}
>0 &\text{ if }y\in(0,m_{p,a}), \\
=0 &\text{ if }y=m_{p,a}, \\
<0 &\text{ if }y\in(m_{p,a},\infty).
\end{cases}
$$
Moreover,
$$(1)\qquad g_{1;p,a}(y):=g'_{p,a}(y)(1 + y)^2 (a + y)^2 y^{1 - p}=(a+y) (p y+p-y)-y (y+1)$$
equals $g'_{p,a}(y)$ in sign. So, if $y=m_{p,a}$, then $g_{1;p,a}(y)=0$. This and $(1)$ imply
$\dfrac{\partial}{\partial a}g_{1;p,a}(y)=p y+p-y=\dfrac{y (y+1)}{a+y}>0$ at $y=m_{p,a}$. It follows that the mode $m_{p,a}$ is increasing in $a$.
Take now any $a_1$ and $a_2$ such that $0\le a_1<a_2\le1$. Then $m_{p,a_1}<m_{p,a_2}$, by the monotonicity of $m_{p,a}$ in $a$. Moreover, by $(0)$, the condition
$$(2)\qquad m_{p,a_1}<y<m_{p,a_2}$$
is equivalent to $g_{1;p,a_1}(y)<0<g_{1;p,a_2}(y)$, which in turn reduces to
$$(3)\qquad\frac{2 y}{1 + y} < p < \frac{1 + 2 y}{1 + y} \quad\&\quad 0 \le a_1 < a_{p, y} \quad \&\quad
a_{p, y} < a_2 \le 1,$$
where $a_{p, y}:=\frac{y - p y + 2 y^2 - p y^2}{p - y + p y}\in(0,1)$, and at that $p - y + p y>0$.
Let now $\mathrm{D}g_{p,a}(y) := g'_{p,a}(y) (1 + y)^2 y^{1 - p}$ and $\mathrm{D}\mathrm{D}g_{p,a}(y) := g''_{p,a}(y) (1 + y)^2 y^{2 - p}$. Then, by Remark 1 in Kemperman-modes, the result in question holds with the exponent $3/2$ replaced by $p$ iff
$$d_{p,a_1,a_2}(y):=
\begin{vmatrix}
\mathrm{D}g_{p,a_1}(y) & \mathrm{D}\mathrm{D}g_{p,a_1}(y) \\
\mathrm{D}g_{p,a_2}(y) & \mathrm{D}\mathrm{D}g_{p,a_2}(y)
\end{vmatrix}
\;\frac{(a_1 + y)^3 (a_2 + y)^3}{y (1 + y)^2 (a_2 - a_1)} $$
$$=\left(a_1+a_2\right) (p-1) y (p y+p-2 y)+a_1 a_2 p (p y+p-y+1)+(p-2) y^2 (p y+p-3 y-1)\ge0$$
whenever $(2)$ or, equivalently, $(3)$ holds.
Since $d_{p,a_1,a_2}(y)$ is affine in $a_1$ and in $a_2$, it suffices to check the latter inequality just for the pairs $(a_1,a_2)\in\{(0,a_{p,y}),(0,1),(a_{p,y},a_{p,y}),(a_{p,y},1)\}$, assuming also $\frac{2 y}{1 + y} < p < \frac{1 + 2 y}{1 + y}$. This checking can be done comparatively quickly (say using the mentioned command Reduce[]), which shows that $d_{p,a_1,a_2}(y)>0$ for $(a_1,a_2)\in\{(0,a_{p,y}),(a_{p,y},a_{p,y}),(a_{p,y},1)\}$, whereas $d_{p,0,1}(y)<0$ for all $p\in(0,1)$ and all $y\in(0,y_p)$ with $y_p:=\dfrac{\left(\sqrt{2}-\sqrt{1-p} \sqrt{2-p}\,\right) \sqrt{1-p}}{(3-p)\sqrt{2-p}}>0$.
Thus indeed, the unimodality holds precisely for $p\in[1,2)$.
