  A time-stepping $L1$ scheme for solving a time fractional Fokker-Planck
equation of order $\alpha \in (0, 1)$, with a general driving force, is
investigated. A stability bound for the semi-discrete solution is obtained for
$\alpha\in(1/2,1)$ {via a novel and concise approach.} Our stability estimate
is $\alpha$-robust in the sense that it remains valid in the limiting case
where $\alpha$ approaches $1$ (when the model reduces to the classical
Fokker-Planck equation), a limit that presents practical importance. Concerning
the error analysis, we obtain an optimal second-order accurate estimate for
$\alpha\in(1/2,1)$. A time-graded mesh is used to compensate for the singular
behavior of the continuous solution near the origin. The $L1$ scheme is
associated with a standard spatial Galerkin finite element discretization to
numerically support our theoretical contributions. We employ the resulting
fully-discrete computable numerical scheme to perform some numerical tests.
These tests suggest that the imposed time-graded meshes assumption could be
further relaxed, and we observe second-order accuracy even for the case
$\alpha\in(0,1/2]$, that is, outside the range covered by the theory.
