**1. Introduction**

Differential and integral operators are the basis of mathematical models, and they are also used as a means of understanding the working principles of natural and artificial systems. Therefore, differential and integral equations are of great importance both theoretically and practically. Such equations have a wide range of applications, including in the physical sciences (such as in physics and engineering) as well as in social science. Systems of differential equations, as differential equations, are often used in issues such as theories of elasticity, dynamics, fluid mechanics, oscillation, and quantum dynamics.

Interest in differential and integral operators has led to the exploration of fractional differential and integral operators by examining these issues further in depth. Owing to a question, the origin of fractional calculus arose in a message from Leibniz to L'Hôpital in 1695. Fractional calculus has received attention in recent years due to its ability to simplify numerous physical, engineering, and economics phenomena such as the fluid dynamic traffic model, damping laws, continuum and statistical mechanics, diffusion processes, solid mechanics, control theory, colored noise, viscoelasticity, electrochemistry, and electromagnetism, among others.

Because a variety of solutions of fractional differential equations (FDEs) cannot be found analytically, numerical and approximate methods are needed. There are a lot of techniques that have been studied by many researchers in solving FDEs and the system of such equations numerically. Some of these techniques are the Adomian decomposition method presented by Song and Wang [1], the collocation method, the improved operational matrix method [2–4], the perturbation iteration method introduced by Senol and Dolapçı [ ¸ 5], the computational matrix method illustrated by Khader et al. [6], the differential transform method demonstrated by Ertürk and Momani [7], the variational iteration method, the Laplace transform method given by Gupta et al. [8], and the fractional complex transform method studied by Ghazanfari and Ghazanfari [9], among others. Kilbas et al. [10] inclusively examined fractional differential and fractional integro-differential equations. In addition, numerical solutions of FDEs and the system of such equations have been presented using the Legendre polynomial operational matrix method [11], Bernstein operational matrix method [12], Genocchi operational matrix method [13], Jacobi operational matrix method [14], Chebyshev wavelet operational matrix method [15], polynomial least squares method (PLSM) [16], Legendre wavelet-like operational matrix method (LWPT) [17], and the Genocchi wavelet-like operational matrix method [18].

This paper focuses on the numerical analysis of a system of fractional order differential equations using the Legendre wavelet operational matrix method. The most important advantage of the proposed method is that it presents an understandable procedure to reduce FDEs and the system of such equations to a system of algebraic equations. First, we begin by presenting some basic definitions and fundamental relations in Sections 2 and 3, respectively. Then, in Section 4, the operational matrix of the fractional derivate is natively formulated to linear and nonlinear systems of fractional differential equations. Section 5 presents five illustrative examples that were tested with the introduced method. Finally, the last section includes the conclusions.
