We develop a numerical method for the Levy-Fokker-Planck equation with the fractional diffusive scaling. There are two main challenges. One comes from a two-fold non locality, that is, the need to apply the fractional Laplacian operator to a power law decay distribution. The other comes from long-time/small mean-free-path scaling, which calls for a uniform stable solver. To resolve the first difficulty, we use a change of variable to convert the unbounded domain into a bounded one and then apply Chebyshev polynomial 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 fractional diffusion limit analytically.[-]