Abstract
This manuscript presents a new approximation method for fractional-order Fokker-Planck equations based on Touchard polynomial approximation. We provide new Caputo and extra Caputo pseudo-operational matrices for these polynomials. Then, utilizing mentioned pseudo-operational matrices and an optimal method, the considered equation leads to a system of algebraic equations which can be solved by mathematical software. Finally, we illustrate the advantages of the suggested technique through several numerical examples.