$\newcommand{\Z}{\mathbb Z}$Of course, one can write an expression for the probability in question as a certain sum over the set $\{-1,0,1\}^N$. Given the three parameters $N,M,z$, a very simple expression does not seem likely to exist.
However, looking at the distribution of $X_i-Y_i$, one may note that the probability in question, say $p_{z,N}=p_{z,n,M}$, is the solution to the following boundary value problem for a partial difference equation:
$$p_{y,n}=\tfrac14\,p_{y+1,n-1}+\tfrac14\,p_{y-1,n-1}+\tfrac12\,p_{y,n-1} \tag{1}\label{1}$$
for $(y,n)\in\{1-M,\dots,M-1\}\times\{1,\dots,N\}$,
$$p_{y,0}=1(y<0)\quad \text{ for $y\in\{-M,\dots,M\}$,} \tag{2}\label{2}$$
$$p_{M,n}=0,\quad p_{-M,n}=1 \quad \text{ for $n\in\{1,\dots,N\}$. }\tag{3}\label{3}$$
This makes it easy to compute the solution for any given $N,M,z$.
Mathematica's command RSolve cannot do anything with partial difference equation \eqref{1}. However, clearly, the formula $p_{y,n}=s^y t_s^n$ defines a solution to \eqref{1} for any complex $s\ne0$ and $t_s:=\frac s4+\frac1{4s}+\frac12$. So, if (say) for some natural $k$ and some complex $c_1,\dots,c_k,s_1,\dots,s_k$ the formula
$$p_{y,n}=\sum_{j=1}^k c_j s_j^y t_{s_j}^n \tag{4}\label{4}$$
defines a solution to the boundary conditions \eqref{2} and \eqref{3}, then \eqref{4} actually defines the unique solution to the boundary value problem \eqref{1}--\eqref{3}. More generally, one may want to try to find the solution to \eqref{1}--\eqref{3} in the form
$$p_{y,n}=\int_U \mu(du) s_u^y t_{s_u}^n \tag{4a}\label{4a}$$
for some complex-valued measure $\mu$ on some measurable space $(U,\mathcal U)$ and some measurable function $U\ni u\mapsto s_u\in\Bbb C$.
Another expression for $p_{y,n}$ can be obtained using the reflection principle and the multinomial (here, more specifically, trinomial) formula.
Indeed, for any integer $z\in(-M,M)$,
\begin{equation*}
p_{z,N}=\sum_{n=0}^{N-1}P(|z+S_j|<M\ \forall j\in[n],\ z+S_n=1-M,U_{n+1}=-1),
\end{equation*}
where $S_n:=\sum_{i=1}^n U_i$ (so that $S_0=0$), $U_i:=X_i-Y_i$, and $[n]:=\{1,\dots,n\}$. The $U_i$'s are i.i.d., with $P(U_i=\pm1)=1/4$ and $P(U_i=0)=1/2$.
So,
\begin{equation*}
p_{z,N}=\frac14\sum_{n=0}^{N-1}q_n, \tag{10}\label{10}
\end{equation*}
where
\begin{equation*}
q_n:=q_{M,z,n}:=P(S_j\in(a_-,a_+)\ \forall j\in[n],\ S_n=b),
\end{equation*}
where
\begin{equation*}
a_\pm:=\pm M-z,\quad b:=1-M-z.
\end{equation*}
In turn, by multiple reflection of the symmetric random walk $(S_j)$ with jumps in the set $\{-1,1\}$ (cf. e.g. Proposition 4),
\begin{equation*}
q_n=\sum_{k\in\Z}(-1)^kP(S_n=b_k), \tag{20}\label{20}
\end{equation*}
where
\begin{equation*}
b_k:=(-1)^k b+\frac{1-(-1)^k}2(a_++a_-)+k(a_+-a_-) \\
=(2k-(-1)^k)M-z+(-1)^k;
\end{equation*}
note that only finitely many summands in \eqref{20} are nonzero -- namely, the summands with $k$ such that $|b_k|\le n$.
Further, $P(S_n=u)=0$ is $|u|>n$, and for integers $u\in[-n,n]$, by the trinomial formula,
\begin{align*}
P(S_n=u)&=\sum_{(n_0,n_1,n_{-1})\in T_{n,u}}\frac{n!}{n_0! n_1! n_{-1}!}
\Big(\frac12\Big)^{n_0}\Big(\frac14\Big)^{n-n_0} \\
&=2^{-n-|u|}n!\sum_{j=0}^{j_{n,u}}\frac{2^{-2j}}{(n-|u|-2j)!(|u|+j)! j!}, \tag{30}\label{30}
\end{align*}
where $j_{n,u}:=\lfloor\frac{n-|u|}2\rfloor$ and
$T_{n,u}$ is the set of all triples $(n_0,n_1,n_{-1})$ of nonnegative integers such that $n_0+n_1+n_{-1}=n$ and $n_1-n_{-1}=u$.
Collecting \eqref{10}, \eqref{20}, and \eqref{30}, we express $p_{z,N}$ as a triple sum, in $n,k,j$;
\begin{equation}
p_{z,N}=\frac14\sum_{n=0}^{N-1} 2^{-n}n!
\sum_{k\in\Z}(-1)^k
2^{-|b_k|}\sum_{j=0}^{j_{n,b_k}}\frac{2^{-2j}}{(n-|b_k|-2j)!(|b_k|+j)! j!};
\end{equation}
note that, if $|b_k|>n$, then $j_{n,b_k}<0$ and hence, by the standard convention, $\sum_{j=0}^{j_{n,b_k}}\ldots=0$.
