I am currently looking for a discrete version of the parabolic Harnack inequality in the following "$L^1$ to $L^\infty$" form:

If $u(t,x)\geq 0$ is a (say, smooth) subsolution of \begin{equation} \partial_t u -\Delta u \leq C u \tag{1} \end{equation} for some constant $C>0$ then \begin{equation} u(t_0,x_0)\leq \frac{C'}{ R^{d+2}}\int _{t_0-R^2}^{t_0}\int_{B_R(x_0)}u(t,x)dt\,dx \tag{2} \end{equation} Here $d>0$ is the space dimension so the $R^{d+2}$ factor comes from the correct parabolic scaling. Note that I'm only interested in a local version on balls, I don't care so much about boundary conditions.

I'm trying to figure out how to prove an equivalent statement for any numerical scheme. I'm not sure yet which scheme exactly I'm going to need so allow me to remain vague at this stage, right now I'm looking for a reference to get started. To fix ideas, assume that I'm using one-dimensional finite differences with mesh sizes $\Delta t,\Delta x\ll 1$, write $x_i=i\Delta t$ for $i\in\mathbb Z$ and $t^n=n\Delta t$ for $n\in\mathbb N$. Then, for the approximation $$ u^n_i\approx u(t^n,x_i)\qquad n\in\mathbb N,\,i\in\mathbb Z, $$ a natural (explicit in time) discrete version of our differential inequality (1) takes the form $$ \frac{u^{n+1}_i-u_i^n}{\Delta t}-\frac{u^n_{i+1}-2u^n_{i}+u^n_{i-1}}{\Delta x^2}\leq C u_i^n. $$

Question:can we prove a discrete version of (2) in this simple case? (or for any variant thereof) Is anyone aware of any result along these lines to get me started?

Any help will be greatly appreciated!