  In this paper, we develop a numerical method for the L\'evy-Fokker-Planck
equation with the fractional diffusive scaling. There are two main challenges.
One comes from a two-fold nonlocality, that is, the need to apply the
fractional Laplacian operator to a power law decay distribution. The other
arises from long-time/small mean-free-path scaling, which introduces stiffness
to the equation. To resolve the first difficulty, we use a change of variable
to convert the unbounded domain into a bounded one and then apply the Chebyshev
polynomial based pseudo-spectral method. To treat the multiple scales, we
propose an asymptotic preserving scheme based on a novel micro-macro
decomposition that uses the structure of the test function in proving the
fractional diffusion limit analytically. Finally, the efficiency and accuracy
of our scheme are illustrated by a suite of numerical examples.
