  In this paper, we introduce and analyse numerical schemes for the homogeneous
and the kinetic L\'evy-Fokker-Planck equation. The discretizations are designed
to preserve the main features of the continuous model such as conservation of
mass, heavy-tailed equilibrium and (hypo)coercivity properties. We perform a
thorough analysis of the numerical scheme and show exponential stability and
convergence of the scheme. Along the way, we introduce new tools of discrete
functional analysis, such as discrete nonlocal Poincar\'e and interpolation
inequalities adapted to fractional diffusion. Our theoretical findings are
illustrated and complemented with numerical simulations.
