In this paper, we introduce and analyse numerical schemes for the homogeneous and the kinetic Lévy-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. Along the way, we introduce new tools of discrete functional analysis, such as discrete nonlocal Poincaré and interpolation inequalities adapted to fractional diffusion. Our theoretical findings are illustrated and complemented with numerical simulations.