  We develop a unified and easy to use framework to study robust fully discrete
numerical methods for nonlinear degenerate diffusion equations $$ \partial_t
u-\mathfrak{L}^{\sigma,\mu}[\varphi(u)]=f \quad\quad\text{in}\quad\quad
\mathbb{R}^N\times(0,T), $$ where $\mathfrak{L}^{\sigma,\mu}$ is a general
symmetric diffusion operator of L\'evy type and $\varphi$ is merely continuous
and non-decreasing. We then use this theory to prove convergence for many
different numerical schemes. In the nonlocal case most of the results are
completely new. Our theory covers strongly degenerate Stefan problems, the full
range of porous medium equations, and for the first time for nonlocal problems,
also fast diffusion equations. Examples of diffusion operators
$\mathfrak{L}^{\sigma,\mu}$ are the (fractional) Laplacians $\Delta$ and
$-(-\Delta)^{\frac\alpha2}$ for $\alpha\in(0,2)$, discrete operators, and
combinations. The observation that monotone finite difference operators are
nonlocal L\'evy operators, allows us to give a unified and compact {\em
nonlocal} theory for both local and nonlocal, linear and nonlinear diffusion
equations. The theory includes stability, compactness, and convergence of the
methods under minimal assumptions -- including assumptions that lead to very
irregular solutions. As a byproduct, we prove the new and general existence
result announced in \cite{DTEnJa17b}. We also present some numerical tests, but
extensive testing is deferred to the companion paper \cite{DTEnJa18b} along
with a more detailed discussion of the numerical methods included in our
theory.
