Anomalous diffusion is one of the most ubiquitous phenomena in nature, and it is present in a wide variety of physical situations, for instance, transport of fluid in porous media, diffusion of plasma, diffusion at liquid surfaces, etc. The fractional approach proved to be highly effective in a rich variety of scenarios such as continuous time random walk models, generalized Langevin equations, or the generalized master equation. To investigate the subdiffusion of anomalous diffusion, it would be useful to study a time fractional Fokker–Planck equation. In this paper, firstly the time fractional, the sense of Riemann–Liouville derivative, Fokker–Planck equation is transformed into a time fractional ordinary differential equation (FODE) in the sense of Caputo derivative by discretizing the spatial derivatives and using the properties of Riemann–Liouville derivative and Caputo derivative. Then combining the predictor–corrector approach with the method of lines, the algorithm is designed for numerically solving FODE with the numerical error O(k^min{1+2α,2})+O(h²), and the corresponding stability condition is got. The effectiveness of this numerical algorithm is evaluated by comparing its numerical results for α=1.0 with the ones of directly discretizing classical Fokker–Planck equation, some numerical results for time fractional Fokker–Planck equation with several different fractional orders are demonstrated and compared with each other, moreover for α=0.8 the convergent order in space is confirmed and the numerical results with different time step sizes are shown.