We develop a numerical method for the Levy-Fokker-Planckequation with the fractional diffusive scaling. There are two main challenges. One comes from a two-fold nonlocality, that is, the need to apply the frac-tional Laplacian operator to a power law decay distribution. The other comes from long-time/small mean-free-path scaling, which calls for a uniform stablesolver. To resolve the first difficulty, we use a change of variable to convert theunbounded domain into a bounded one and then apply Chebyshev polyno-mial based pseudo-spectral method. To resolve the second issue, we propose an asymptotic preserving scheme based on a novel micro-macro decomposition that uses the structure of the test function in proving the fractionaldiffusion limit analytically.