The goal of this work is to develop a high-order numerical technique to solve a two-dimensional fractional neutron diffusion (FND) model describing dynamical behavior of a lead-cooled fast reactor. This method is based on L1 technique for temporal discretization and a high-order compact Alternating Direction Implicit (ADI) finite difference scheme for the spatial discretization. Numerical simulations are conducted to demonstrate the applicability of the method. We calculate the neutronic flux in the core and examine how the order of the temporal-fractional derivative influences the neutronic flux. The results obtained with FND model for different values of anomalous diffusion coefficient are compared against those obtained with P1 model and classical neutron diffusion equation. In addition, the computational time (CPU time) is provided in order to justify the computational efficiency of the proposed method.