**6. Conclusions**

A practical matrix approach based on novel (orthogonal) Bessel polynomials is presented to solve multi-order fractional-order differential equations (MOFDEs). Using the matrix representations of the generalized Bessel polynomials and their derivatives with the aid of collocation points, the scheme transforms MOFDEs to a fundamental matrix equation, which corresponds to a system of (non)linear algebraic equations. To assess the efficiency and accuracy of the presented technique, several numerical examples with initial and boundary conditions are investigated. Comparisons with the exact solutions and with various alternative numerical simulations and experimental measurements have also been made. Based on the experiments, it is found that the numerical approximations are in an excellent agreement, which demonstrate the reliability and the grea<sup>t</sup> potential of the presented technique.

**Author Contributions:** Conceptualization, M.I. and C.C.; methodology, M.I. and C.C.; software, M.I.; validation, M.I. and C.C.; formal analysis, M.I. and C.C.; investigation, M.I. and C.C.; writing—original draft preparation, M.I. ; writing—review and editing, M.I. and C.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
