If $b$ is fixed (which is what was actually needed in the linked problem at incr-decr parametric integral), by rescaling $z:=y/b$ without loss of generality (wlog) $b=1$.
The crucial ingredient in the proof is Theorem 4 in Kemperman-modes, which 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
$$(**)\qquad 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 for a fixed $b$.
In contrast, the case with $b$ allowed to freely vary in $(0,\infty)$ has just been answered, negatively, by Douglas Zare.
Addendum: I got curious as to how essential the exponent $3/2$ in the definition $(*)$ in question is. It turns out that the unimodality result (for a fixed $b$) 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)$.
