In this study, the Caputo fractional derivative is used to define time fractional Fokker-Planck equation in an unbounded domain. To solve this equation, the Jacobi polynomials together with the tanh-Jacobi functions are employed. The operational matrices of the classical and fractional derivatives of these basis functions are obtained to use them in constructing a numerical method for the expressed equation. In the proposed method, the introduced basis functions are used simultaneously to approximate the equation’s unknown solution. More precisely, the shifted Jacobi polynomials are applied to approximate the solution in the temporal direction and the tanh-Jacobi functions are utilized to approximate the solution in the spatial direction. By substitute the expressed approximation into the equation and employing the introduced operational matrix, solving the problem under consideration transforms into solving an algebraic system of equations, which can be solved easily. The accuracy and efficiency of the presented method are investigated numerically by solving some numerical examples. The reported results confirms the high accuracy of the established method.