  The fractional Fokker-Planck system with multiple internal states is derived
in [Xu and Deng, Math. Model. Nat. Phenom., $\mathbf{13}$, 10 (2018)], where
the space derivative is Laplace operator. If the jump length distribution of
the particles is power law instead of Gaussian, the space derivative should be
replaced with fractional Laplacian. This paper focuses on solving the two state
Fokker-Planck system with fractional Laplacian. We first provide a priori
estimate for this system under different regularity assumptions on the initial
data. Then we use $L_1$ scheme to discretize the time fractional derivatives
and finite element method to approximate the fractional Laplacian operators.
Furthermore, we give the error estimates for the space semidiscrete and fully
discrete schemes without any assumption on regularity of solutions. Finally,
the effectiveness of the designed scheme is verified by numerical experiments.
