$\def\RR{\mathbb{R}}$This isn't the right answer for your class, because it works for a level of generality you don't need, but last summer I found what I think is the right way to prove this result for a general bipartite planar graph or, working harder, for any planar graph and I thought I'd record it here. I'm sure this isn't original, but I didn't find it in a quick literature search.
Let $G$ be a bipartite graph. We'll define an immersion of $G$ to $\RR^2$ to be a map $\phi: G \to \RR^2$ which is injective except that $\phi$ is allowed to make edges cross transversely in their interiors; $\phi$ is not permitted to create triple crossings. For example, we can send the vertices of $G$ to generic points of $\RR^2$ and join them by line segments.
Assign a weight $x_e$ to each edge of $G$. For a perfect matching $M$ of $G$, define $\mathrm{cross}(M, \phi)$ to be the number of places where edges of $M$ cross. Define the matching polynomial of $(G, \phi)$ to be
$$\mu(G, \phi) = \sum_M (-1)^{\mathrm{cross}(M, \phi)} \prod_{e \in M} x_e$$
where the sum is over all perfect matchings of $G$.
We claim that there is a way to take the adjacency matrix of $G$ and replace the entry corresponding to $e$ by $\pm x_e$ to make a matrix $A$ with $\mu(G, \phi) = \det A$. In particular, if $G$ is a planar graph, we can get the generating function for matchings with no signs.
Proof First, suppose that all of the black vertices lie on the line $y=0$ and all the white vertices lie on the line $y=1$, in the same order as the rows and columns of the matrix, and joined by line segments. Then we can put $x_e$ in all positions of $A$, and $(-1)^{\mathrm{cross}(\phi, M)}$ is the sign of the permutation $M$.
Now, consider varying $\phi$ while preserving $G$ and we claim that the result stays true. If we vary $\phi$ generically, we can get from any immersion to any other while passing through singularities of the following sorts:
A triple crossing of edges. The same signs work on both sides of the triple crossing.
Two edges become tangent. On one side of the singularity, they cross $m$ times, on the other side they cross $m+2$ times. The same signs work on both sides of the tangency.
An edge $e$ forms a cusp. On one side of the cusp, it passes through itself, on the other side it doesn't. Negate the sign of $x_e$.
$\phi$ sends a vertex into the interior of an edge $e$. Negate the sign of $x_e$.
So in every case, as the topology of $\phi$ changes, we can change the signs in the matrix. $\square$
The analogous argument works with general graphs and Pfaffians, starting with the vertices all lying on a circle.
Your difference operator $\Delta^\mu$ has the property that it commutes with the shift
operator $Ef(x)=f(x+1)$. It therefore falls in the purview of the work of Rota, et al., on
finite operator calculus. In particular, Theorem 2 (First Expansion Theorem) of G.-C. Rota, D. Kahaner, and A. Odlyzko, On the foundations of combinatorial theory. VIII. Finite operator calculus, J. Math. Anal. Appl. 42 (1973), 684-760, gives the discrete Taylor expansion you are looking for. For reasons of formal convergence, you need the additional condition that $\Delta^\mu(x)$ is a nonzero constant.
Best Answer
Have you considered using the matrix-tree theorem [2, 3] which counts spanning trees instead of domino tilings?
For a rectangular graph (or any product of intervals) we can engineer a graph whose eigenvalues multiply to $K_r$ (not that I have worked it out).
There is a bijection by Temperley which connects domino tilings and spanning trees (since both involve planar graphs). I don't know of higher dimensional Temperley bijection.
Tzeng + Wu [Spanning Trees on Hypercubic Lattices and Non-Orientable Surfaces] (also works out Möbius strip and Klein bottle).
For a modern discussion of Kasteleyn's theorem look at
Variations on a theme of Kasteleyn, with application to the totally nonnegative Grassmannian