This research presents a numerical approach to obtain the approximate solution of the n-dimensional cohomological equations of fractional order in continuous-time dynamical systems. For this purpose, the $ n $-dimensional fractional M\”{u}ntz-Legendre polynomials (or n-DFMLPs) are introduced. The operational matrix of the fractional Riemann-Liouville derivative is constructed by employing n-DFMLPs. Our method transforms the cohomological equation of fractional order into a system of algebraic equations. Therefore, the solution of that system of algebraic equations is the solution of the associated cohomological equation. The error bound and convergence analysis of the applied method under the $ L^{2} $-norm is discussed. Some examples are considered and discussed to confirm the efficiency and accuracy of our method.