A novel method is presented for the simulation of a discrete state space, continuous time Markov process subject to fractional diffusion. The method is based on Lie-Trotter operator splitting of the diffusion and reaction terms in the master equation. The diffusion term follows a multinomial distribution governed by a kernel that is the discretized solution of the fractional diffusion equation. The algorithm is validated and simulations are provided for the Fisher-KPP wavefront. It is shown that the wave speed is dictated by the order of the fractional derivative, where lower values result in a faster wave than in the case of classical diffusion. Since many physical processes deviate from classical diffusion, fractional diffusion methods are necessary for accurate simulations.
