In this research, a spectral method is considered for solving d –dimensional nonlinear time–fractional Fokker–Planck equations in an unbounded domain and enhanced for solving the fractional Black-Scholes equations. First, separately for \(d=1,2\) , and 3, this equation is spatially discretized using the Hermite–Galerkin spectral method, and a system of ordinary time-fractional differential equations (TFDE) is obtained. Then, after applying an eigenvalue decomposition technique, a set of independent ordinary TFDEs is obtained, which can be temporally discretized using a Petrov–Galerkin method based on generalized Jacobi functions or solved analytically. Furthermore, a novel approach which is referred to as the extended eigenvalue decomposition and requires less memory and computational complexity compared to other existing approaches, is introduced for discretizing the equation in any dimension. While most of the existing computational methods solve the equations on a truncated domain, which can be ineffective for equations whose solutions decay slowly, we solve the equations on an unbounded domain directly. It is completely natural to choose Hermite functions as basis functions for the models defined on the whole real line, because classical mesh-based methods are ineffective in an unbounded domain. Also, using Hermite functions leads to sparse matrices, and for equations without correlation terms leads to symmetric matrices. Since the exact solutions of fractional differential equations may exhibit a weak singularity at \(t=0\) , spectral methods based on generalized Jacobi functions can significantly improve the accuracy. Also, the convergence criteria of the method for solving high-dimensional nonlinear models and the error estimate for the method are investigated. Simple implementation, high accuracy, flexibility to higher dimensions, applicability to models in finance and physics, and the ability to deal with non-homogeneous boundaries at infinity are advantages of this approach.