Consider the standard 5 point different approximation (centered difference for both the gradient and divergence operators) for the variable coefficient Poisson equation $$-\nabla \cdot(a\nabla v)=f$$ with Dirichlet boundary conditions, in a 2 dimensional rectangular region. We assume that $a(x,y)\geq a_0 >0$. The approximate solution $\{u_{i,j}\}$ satisfies a linear system $Au=b$.
Question: State and prove the maximum principle for the numerical solution $u_{i,j}$.
My attempt: I am a bit unsure if the minus sign changes this problem. Maybe the $\geq$ should be a $\leq$ or something? Any help on this would be appreciated. I am also not sure if I am incorporating the variable coefficients correctly or not. Do I need to include them as I did or do I only consider the $u_{i,j}$?
Statement: Assume that $-\nabla \cdot(a\nabla u_{i,j}) \geq 0$ for all $i=1,2,...,m$ and for all $j=1,2,...,n$. Let M be the maximum in the rectangular grid. Then, if $a_{i0,j0}u_{i0,j0}=M$ on any interior grid point, then $u_{i,j}=M$ for all interior and boundary points.
Proof: We have that $$M=-\nabla \cdot(a_{i0,j0}\nabla u_{i0,j0})\geq 0$$ for some point on the interior of the grid. Let h be delta x and k be delta y. $$\Rightarrow 0 \leq \frac{1}{h^2}\bigg(a_{i0-1,j0}u_{i0-1,j0}-2a_{i0,j0}u_{i0,j0} +a_{i0+1,j0}u_{i0+1,j0}\bigg) + \frac{1}{k^2}\bigg(a_{i0,j0-1}u_{i0,j0-1}-2a_{i0,j0}u_{i0,j0} +a_{i0,j0+1}u_{i0,j0+1}\bigg)$$
$$\Rightarrow 0 \leq \frac{1}{h^2}\bigg(a_{i0-1,j0}u_{i0-1,j0}-2M +a_{i0+1,j0}u_{i0+1,j0}\bigg) + \frac{1}{k^2}\bigg(a_{i0,j0-1}u_{i0,j0-1}-2M +a_{i0,j0+1}u_{i0,j0+1}\bigg)$$
So we must have the neighboring terms also equal to M or else the RHS will be smaller than 0. Repeat the argument for the neighbors until the entire region is shown to be M.