Summary: We develop a finite volume method to numerically solve the \(N\)-dimensional time fractional Fokker-Planck equation \[ \frac{\partial ^\alpha \omega }{\partial t^\alpha }=k_\alpha \Delta \omega -\mathop{\sum}\limits _{k=1}^N\frac{\partial (f^{(k)}\omega )}{\partial x_k}, \] where \(\frac{\partial ^\alpha \omega }{\partial t^\alpha }\) is the Caputo fractional derivative of order \(\alpha\) with \(0<\alpha <1\). We theoretically prove that the method is unconditional stable and the error of the numerical solution is \(O(h+\Delta t^{2-\alpha })\) which can be improved to \(O(h^2+\Delta t^{2-\alpha })\) when the space grid is sufficiently fine, where \(h\) and \(\Delta t\) denote space and time grid sizes, respectively. Numerical tests are conducted to support our theoretical analysis.