We assume all solutions to be bounded here!
Let $y_{+},y_{-} \in C_c^{\infty}$ be two positive functions.
If we then consider the heat equation
$$\partial_t u(t,x) = \Delta u(t,x)$$ for two different initial conditions
$u_1(0,x) = y_+$ and $u_2(0,x)=y_+-y_-$
then it is clear that $u_1(t,x) \ge u_2(t,x)$ for all $t,x.$
Now, let us consider the nonlinear heat equation
$$\partial_t u(t,x) = \Delta u(t,x)-\vert u(t,x) \vert^2$$
and assume that there exist two smooth solutions $u_1,u_2$ with initial data
$u_1(0,x) = y_+$ and $u_2(0,x)=y_+-y_-$
respectively.
Does this imply $u_1(t,x) \ge u_2(t,x)$?
This sounds plausible, as positive initial data stay positive and negative ones stay negative under evolution. Moreover, at the interface of positive and negative parts, the nonlinearity is small. Still, the nonlinearity makes it hard to say how this linear combination of initial data gets mixed up in the time evolution.
My thoughts
By Duhamel's formula
$$u_1(t,x) = e^{t\Delta}u_1(x) - \int_0^t e^{(t-s) \Delta}\vert u_1(s,x) \vert^2 \ ds.$$
$$u_2(t,x) = e^{t\Delta}u_2(x) - \int_0^t e^{(t-s) \Delta}\vert u_2(s,x) \vert^2 \ ds.$$
Subtracting the two from each other yields
$$u_1(t,x)-u_2(t,x) = e^{t\Delta}y_{-}(x) -\int_0^t e^{(t-s) \Delta}\left(\vert u_1(s,x) \vert^2-\vert u_2(s,x) \vert^2\right) \ ds.$$
Indeed the term $e^{t\Delta}y_{-}(x)$ is positive. Not sure how to conclude anything from the second term though. But it seems like for small $t$ this term is close to zero.
Observation
If instead we take initial conditions $u_1(0,x) = y(x)$ and $u_2(0,x)= \widetilde{y}(x)$ with $0\le \widetilde{y} \le y$ we seem to get $u_1(t,x) \ge u_2(t,x)$ by the above formula.