In this paper, radial basis functions (RBFs) method is proposed for numerical solution of the Liouville–Caputo time- and Riesz space-fractional Fokker–Planck equation with a nonlinear source term. The left-sided and the right-sided Riemann–Liouville fractional derivatives of RBFs are computed and utilized to approximate the spatial fractional derivatives of the unknown function. Also, the time-fractional derivative is discretized by the high order formulas introduced in [J. X. Cao, C. P. Li and Y. Q. Chen, High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (II), Fract. Calculus Appl. Anal. 18(3) (2015) 735–761]. In each time step, via a collocation method, the computations of fractional Fokker–Planck equation are reduced to a linear or nonlinear system of algebraic equations. Several numerical examples are included to demonstrate the applicability, accuracy and stability of the method. Some comparisons are made with the existing results.