  We investigate quasilinear discrete PDEs $\partial_t u = \Delta^N \varphi(u)+
Kf(u)$ of reaction-diffusion type with nonlinear diffusion term defined on an
$n$-dimensional unit torus discretized with mesh size $\tfrac1N$ for $N\in
{\mathbb N}$, where $\Delta^N$ is the discrete Laplacian, $\varphi$ is a
strictly increasing $C^5$ function and $f$ is a $C^1$ function. We establish
$L^\infty$ bounds and space-time H\"older estimates, both uniform in $N$, of
the first and second spatial discrete derivatives of the solutions. In the
equation, $K>0$ is a large constant and we show how these estimates depend on
$K$. The motivation for this work stems originally from the study of
hydrodynamic scaling limits of interacting particle systems.
  Our method is a two steps approach in terms of the H\"older estimate and
Schauder estimate, which is known for continuous parabolic PDEs. We first show
the discrete H\"older estimate uniform in $N$ for the solutions of the
associated linear discrete PDEs with continuous coefficients, based on the Nash
estimate. We next establish the discrete Schauder estimate for linear discrete
PDEs with uniform H\"older coefficients. The link between discrete and
continuous settings is given by the polylinear interpolations. Since this
operation has a non-local nature, the method requires proper modifications.
  We also discuss another method based on the study of the corresponding
fundamental solutions.
