Let
$$
\mathcal{L}\equiv\sum_{i,j}^{n}a_{i,j}\frac{\partial^{2}}{\partial x_{i}\,\partial x_{j}}+\sum_{i}^{n}b_{i}\frac{\partial}{\partial x_{i}}+c
$$
be uniformly elliptic on $\Omega\equiv\mathbb{R}^{n}\times\left(0,T\right)$.
Let
$$
\mathcal{M\equiv L}-\frac{\partial}{\partial t}
$$
Consider the parabolic equality
$$
\mathcal{M}u=f\text{ on }\Omega.
$$
Lastly, suppose $u\left(\mathbf{x},0\right)$ is monotone increasing
in $\mathbf{x}$ and that $f\left(\mathbf{x},t\right)$ is monotone
decreasing in $\mathbf{x}$ for all times $t$.
A function $g:\mathbb{R}^n \rightarrow \mathbb{R}$ is monotone increasing if for all $\mathbf{x},\mathbf{y}\in\mathbb{R}^n$, $g\left(\mathbf{x}\right)\leq g\left(\mathbf{y}\right)$ whenever $\mathbf{x}\leq \mathbf{y}$. Monotone decreasing is defined similarly.

First, consider the case of $a_{i,j}$, $b_{i}$, $c$ being functions of $t$ only. Define $v\left(\mathbf{x},t\right)=u\left(\mathbf{x}-\epsilon,t\right)$ (where $\epsilon \in \mathbb{R}^n$ with $\epsilon > 0$)
and note that
$$
\mathcal{M}v=f\left(\mathbf{x}-\epsilon,t\right)\text{ on }\Omega.
$$
Let $w=u-v$ so that 
$$
\mathcal{M}w=f\left(\mathbf{x},t\right)-f\left(\mathbf{x}-\epsilon,t\right)\leq0\text{ on }\Omega.
$$
Since $u$ is monotone at time zero, we have $$w\left(\mathbf{x},0\right)=u\left(\mathbf{x},0\right)-v\left(\mathbf{x},0\right)=u\left(\mathbf{x},0\right)-u\left(\mathbf{x}-\epsilon,0\right)\geq0 \text{ on } \mathbb{R}^n.$$
These last  two inequalities yield $w\geq0$ (and hence $u$ is monotone in $\mathbf{x}$) everywhere in $\Omega$.

However, this argument is rather restrictive, as it requires the assumption
that the coefficients $a_{i,j}$, $b_{i}$, and $c$ are constant in space (i.e. w.r.t. $\mathbf{x}$). <strike>Does anyone know of any relaxations to the above?</strike> (see edit below)

**EDIT:** Some things were forgotten that are important: We need to assume the boundedness of the coefficients in the operator, along with $$\liminf_{\left|\mathbf{x}\right|\rightarrow\infty} \mathbf{u} \geq 0$$ for all times.

**EDIT:** I have found a potentially interesting relaxation. Suppose $\mathbf{u}$ is convex in $\mathbf{x}$ on $\Omega$. Let $A\in\mathbb{R}^{n\times n}$
be the matrix with $a_{j,k}$ in the $j^{\text{th}}$ row, $k^{\text{th}}$
column. Since $\mathcal{L}$ is assumed to be uniformly elliptic,
$A$ is positive semidefinite. Denote this by $A\succeq0$. Now, relax the requirement on the coefficients $a_{j,k}$ as follows:
for each $i\in\mathcal{S}$, $A\left(\mathbf{x}+\epsilon,t\right)-A\left(\mathbf{x},t\right)\succeq0$
on $\Omega$ for all $\epsilon>0$. Note that in the one-dimensional
case (i.e. $n=1$), this requirement becomes $a\equiv a_{1,1}$ is monotone
increasing in space. Noting that
$$
\left(A\left(\mathbf{x}+\epsilon,t\right)-A\left(\mathbf{x},t\right)\right)\circ\left[\nabla_{\mathbf{x}}u_{i}\right]\left(\mathbf{x},t\right)\succeq0
$$
where $\nabla$ denotes the Hessian operator and $\circ$ the Hadamard (entrywise) product, we get
$$
\sum_{i,j}^{n}a_{i,j}\left(\mathbf{x},t\right)\frac{\partial^{2}u}{\partial x_{i}\,\partial x_{j}}\leq\sum_{i,j}^{n}a_{i,j}\left(\mathbf{x}+\epsilon,t\right)\frac{\partial^{2}u}{\partial x_{i}\,\partial x_{j}}.
$$
Using the above, we arrive once again at $\mathcal{M}w\leq 0$. This brings me to a more interesting question...

**Are there better relaxations than the above when $u\left(\mathbf{x},t\right)$ is convex in $\mathbf{x}$?**