This paper presents a practical numerical method, an implicit finite-difference scheme for solving a two-dimensional time-space fractional Fokker–Planck equation with space–time depending on variable coefficients and source term, which represents a model of a Brownian particle in a periodic potential. The Caputo derivative and the Riemann–Liouville derivative are considered in the temporal and spatial directions, respectively. The Riemann–Liouville derivative is approximated by the standard Grünwald approximation and the shifted Grünwald approximation. The stability and convergence of the numerical scheme are discussed. Finally, we provide a numerical example to test the theoretical analysis. © 2021 by the authors. Licensee MDPI, Basel, Switzerland.