(Look at page p68-69 in Fulton, Harris Representation Theory) for an explanation of what $\epsilon$ is. http://www.math.ucsb.edu/~bigelow/math227b/fultonharris.pdf
Proposition: $p$ is an odd prime. Let $X \in GL_{2}(\mathbb{Z}_{p})$, $X \notin Z(GL_{2}(\mathbb{Z}_{p}))$, $M \in M_{2}(\mathbb{Z}_{p})$. Then there are precisely $p^{2}$ solutions to $MX-XM=0$. Equivalently, in the Lie algebra $M_{2}(\mathbb{Z}_{p})$ with Lie bracket $[\cdot,\cdot]$ the Lie subalgebra $Y$ of solutions $M$ to $[M,X]=0$ has dimension $2$.
Proof: It is routine to check that these statements are equivalent and that $Y$ is a Lie subalgebra of the Lie algebra $M_{2}(\mathbb{Z}_{p})$. Observe also that for $X,M$ as above and $P \in GL_{2}(\mathbb{Z}_{p})$\, $[M,X]=0 \Longleftrightarrow [P^{-1}MP,P^{-1}XP]=0$. So the dimension of $Y$ is invariant under conjugation of $X$, because if $M_{1},M_{2},M_{3},M_{4}$ is a basis of $M_{2}(\mathbb{Z}_{p})$, then so is $P^{-1}M_{1}P,P^{-1}M_{2}P,P^{-1}M_{3}P,P^{-1}M_{4}P$. So it suffices to check that the statement holds for a representative $X$ from each of the non-central conjugacy classes of $GL_{2}(\mathbb{Z}_{p})$. We use the representatives given in the table on page p68-69 in Fulton, Harris Representation Theory, and have 3 cases:
Case 1: $X=b_{x}=\begin{pmatrix} x & 1 \\ 0 & x \end{pmatrix}$, $x \neq 0$. We have a basis of $M_{2}(\mathbb{Z}_{p})$ given by: \begin{align*}A= \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & B=\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & C=\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} & D=\begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}\end{align*} Now:
\begin{align*}[A,X]&= 0 \\ [B,X] & = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\begin{pmatrix} x & 1 \\ 0 & x \end{pmatrix} - \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\begin{pmatrix} x & 1 \\ 0 & x \end{pmatrix} = \begin{pmatrix} 0 & 2 \\ 0 & 0 \end{pmatrix} \neq 0 \\ [C,X]& = \begin{pmatrix} 0 & 1 \\ 0 & 0\end{pmatrix}\begin{pmatrix} x & 1 \\ 0 & x \end{pmatrix} - \begin{pmatrix} x & 1 \\ 0 & x \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 0 & 0\end{pmatrix} = \begin{pmatrix} 0 & x \\ 0 & 0\end{pmatrix}-\begin{pmatrix} 0 & x \\ 0 & 0\end{pmatrix} = 0 \\ [D,X]&=\begin{pmatrix} 0 & 0 \\ 0 & 1\end{pmatrix}\begin{pmatrix} x & 1 \\ 0 & x \end{pmatrix}-\begin{pmatrix} x & 1 \\ 0 & x \end{pmatrix}\begin{pmatrix} 0 & 0 \\ 0 & 1\end{pmatrix} = \begin{pmatrix} -1 & 0 \\ 0 & 1\end{pmatrix} \neq 0\end{align*}
So $A,C \in Y$ and $dim(Y)\geq 2$. Note that $[B,X],[D,X]$ are linearly independent in $M_{2}(\mathbb{Z}_{p})$, so for $M=\alpha A+\beta B+ \gamma C+ \delta D \in M_{2}(\mathbb{Z}_{p})$, $[M,X]=0\Longleftrightarrow \beta=\delta=0$, so $dim(Y)=2$.
Case 2: $X=c_{x,y}=\begin{pmatrix} x & 0 \\ 0 & y \end{pmatrix}$. $[A,X], [B,X]=0$ in this case. We have:
\begin{align*}[C,X]&=\begin{pmatrix} 0 & 1 \\ 0 & 0\end{pmatrix}\begin{pmatrix} x & 0 \\ 0 & y \end{pmatrix}-\begin{pmatrix} x & 0 \\ 0 & y \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 0 & 0\end{pmatrix} = \begin{pmatrix} 0 & y-x \\ 0 & 0\end{pmatrix} \neq 0 \\ [D,X]&= \begin{pmatrix} 0 & 0 \\ 1 & 0\end{pmatrix}\begin{pmatrix} x & 0 \\ 0 & y\end{pmatrix}-\begin{pmatrix} x & 0 \\ 0 & y\end{pmatrix}\begin{pmatrix} 0 & 0 \\ 1 & 0\end{pmatrix} = \begin{pmatrix} 0 & 0 \\ x-y & 0\end{pmatrix} \neq 0\end{align*} So $A,B \in Y$ and $dim(Y)\geq 2$. Note that $[C,X],[D,X]$ are linearly independent in $M_{2}(\mathbb{Z}_{p})$, so for $M=\alpha A+\beta B+ \gamma C+ \delta D \in M_{2}(\mathbb{Z}_{p})$, $[M,X]=0\Longleftrightarrow \gamma=\delta=0$, so $dim(Y)=2$ in this case.
Case 3: $X=d_{x,y}=\begin{pmatrix} x & \epsilon y \\ y & x \end{pmatrix}$ For this case, we use a slightly different basis for $M_{2}(\mathbb{Z}_{p})$, given by:\begin{align*}A= \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & B=\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & C=\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} & D=\begin{pmatrix} 0 & \epsilon \\ 1 & 0 \end{pmatrix}\end{align*} We have:\begin{align*}[A,X]&=0 \\ [B,X]&=\begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}\begin{pmatrix} x & \epsilon y \\ y & x \end{pmatrix}-\begin{pmatrix} x & \epsilon y \\ y & x \end{pmatrix}\begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix} = \begin{pmatrix}0 & 2\epsilon y \\ -2y & 0\end{pmatrix} \neq 0 \\ [C,X]&=\begin{pmatrix}0 & 1 \\ 0 & 0\end{pmatrix}\begin{pmatrix} x & \epsilon y \\ y & x \end{pmatrix}-\begin{pmatrix} x & \epsilon y \\ y & x \end{pmatrix}\begin{pmatrix}0 & 1 \\ 0 & 0\end{pmatrix}=\begin{pmatrix} y & 0 \\ 0 & -y\end{pmatrix} \neq 0 \\ [D,X]&=\begin{pmatrix} 0 & \epsilon \\ 1 & 0 \end{pmatrix}\begin{pmatrix} x & \epsilon y \\ y & x \end{pmatrix}-\begin{pmatrix} x & \epsilon y \\ y & x \end{pmatrix}\begin{pmatrix} 0 & \epsilon \\ 1 & 0 \end{pmatrix}=0\end{align*}So $A,D \in Y$ and $dim(Y)\geq 2$. Note that $[B,X],[C,X]$ are linearly independent in $M_{2}(\mathbb{Z}_{p})$, so for $M=\alpha A+\beta B+ \gamma C+ \delta D \in M_{2}(\mathbb{Z}_{p})$, $[M,X]=0\Longleftrightarrow \beta=\gamma=0$, so $dim(Y)=2$ in this case, and the statement holds in all 3 cases. $\square$