This paper presents a comprehensive study of non-linear systems of multi-order fractional differential equations from aspects of theory and numerical approximation. In this regard, we first establish the well-posedness of the underlying problem by investigations concerning the existence, uniqueness, and influence of perturbed data on the behavior of the solutions as well as smoothness of the solutions under some assumptions on the given data. Next, from the numerical perspective, we develop and analyze a well-conditioned and high-order numerical approach based on an operational spectral Galerkin method. The main advantage of our strategy is that it characterizes the approximate solution via some recurrence formulas, instead of solving a complex non-linear block algebraic system that requires high computational costs. Moreover, the familiar spectral accuracy is justified in a weighted -norm, and some practical test problems are provided to approve efficiency of the proposed method.