**4. Solving Systems of Fractional Order Differential Equations**

In this section, the Legendre wavelet operational matrix method was implemented to obtain the numerical solution of the system of fractional order differential equations. Consider a system of fractional differential equations as follows:

$$\begin{array}{l} D^{\eta\_1} u\_1(t) = \mathcal{U}\_1(t, \mu\_1, \mu\_2, \dots, \mu\_m), \\ D^{\eta\_2} u\_2(t) = \mathcal{U}\_2(t, \mu\_1, \mu\_2, \dots, \mu\_m), \\ \vdots \\ D^{\eta\_n} u\_m(t) = \mathcal{U}\_m(t, \mu\_1, \mu\_2, \dots, \mu\_m), \end{array} \tag{27}$$

where *Ui* is a linear/nonlinear function of *t*, *u*1, *u*2, ... , *um*, *Dη<sup>i</sup>* is the derivative of *ui* with the order of *η<sup>i</sup>* in the Liouville–Caputo sense and *N* − 1 ≤ *η<sup>i</sup>* < *N*, and they are subjected to the initial conditions

$$\begin{array}{ccccccccc}\mu\_{1}(t\_{0})=\mathsf{u}\_{10} & \frac{d\mathsf{u}\_{1}}{dt}(t\_{0})=\mathsf{u}\_{11}, & \frac{d^{2}\mathsf{u}\_{1}}{dt^{2}}(t\_{0})=\mathsf{u}\_{12}, & \dots, & \frac{d^{n-1}\mathsf{u}\_{1}}{dt^{n-1}}(t\_{0})=\mathsf{u}\_{1(n-1)}\\\mu\_{2}(t\_{0})=\mathsf{u}\_{20}, & \frac{d\mathsf{u}\_{2}}{dt}(t\_{0})=\mathsf{u}\_{21}, & \frac{d^{2}\mathsf{u}\_{2}}{dt^{2}}(t\_{0})=\mathsf{u}\_{22}, & \dots, & \frac{d^{n-1}\mathsf{u}\_{2}}{dt^{n-1}}(t\_{0})=\mathsf{u}\_{2(n-1)}\\\vdots & \vdots & \vdots & \vdots & \vdots\\\mu\_{n}(t\_{0})=\mathsf{u}\_{m0}, & \frac{d\mathsf{u}\_{m}}{dt}(t\_{0})=\mathsf{u}\_{m1}, & \frac{d^{2}\mathsf{u}\_{m}}{dt^{2}}(t\_{0})=\mathsf{u}\_{m2}, & \dots, & \frac{d^{n-1}\mathsf{u}\_{m}}{dt^{n-1}}(t\_{0})=\mathsf{u}\_{m(n-1)}\end{array} \tag{28}$$

First of all, approximating *u*1(*t*), *u*2(*t*), ... , *um*(*t*) and *Dη*1*u*1(*t*), *Dη*2*u*2(*t*), ... , *Dη<sup>n</sup> um*(*t*), we obtain

$$\begin{aligned} \mu\_1(t) &\approx \sum\_{n=0}^{2^k - 1} \sum\_{m=0}^M c\_{1n,m} \psi\_{n,m} = \mathbf{C}\_1^{\ \ \mathbf{T}} \boldsymbol{\upvarphi}(t) \\ \mu\_2(t) &\approx \sum\_{n=0}^{2^k - 1} \sum\_{m=0}^M c\_{2n,m} \psi\_{n,m} = \mathbf{C}\_2^{\ \ \mathbf{T}} \boldsymbol{\upvarphi}(t) \\ &\vdots \\ \mu\_m(t) &\approx \sum\_{n=0}^{2^k - 1} \sum\_{m=0}^M c\_{nm,m} \psi\_{n,m} = \mathbf{C}\_m^{\ \ \ \end{aligned} \tag{29}$$

where *Ci*, *i* = 1, 2, ... , *m* is an unknown vector and *ψ*(*t*) is the vector introduced in Equation (8). If we utilize Equation (17), then we have

$$\begin{array}{l} D^{\eta\_1} u\_1(t) \approx \mathbb{C}\_1{}^T D^{(\eta\_1)} \psi(t) \\ D^{\eta\_2} u\_2(t) \approx \mathbb{C}\_2{}^T D^{(\eta\_2)} \psi(t) \\ \vdots \\ D^{\eta\_n} u\_m(t) \approx \mathbb{C}\_m{}^T D^{(\eta\_n)} \psi(t) \end{array} \tag{30}$$

Substituting Equations (29) and (30) into Equation (27), we obtain

$$\begin{cases} R\_1(t) = \mathbb{C}\_1^T D^{(\eta\_1)} \psi(t) - \mathcal{U}\_1(t, \mathbb{C}\_1^T \psi(t), \mathbb{C}\_2^T \psi(t), \dots, \mathbb{C}\_m^T \psi(t)) \\ R\_2(t) = \mathbb{C}\_2^T D^{(\eta\_2)} \psi(t) - \mathcal{U}\_2(t, \mathbb{C}\_1^T \psi(t), \mathbb{C}\_2^T \psi(t), \dots, \mathbb{C}\_m^T \psi(t)) \\ \vdots \\ R\_m(t) = \mathbb{C}\_m^T D^{(\eta\_n)} \psi(t) - \mathcal{U}\_m(t, \mathbb{C}\_1^T \psi(t), \mathbb{C}\_2^T \psi(t), \dots, \mathbb{C}\_m^T \psi(t)) \end{cases} \tag{31}$$

If *Ui* is a linear function of *<sup>t</sup>*, *<sup>u</sup>*1, *<sup>u</sup>*2, ... , *um*, then we produce 2*k*(*<sup>M</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> *mn* linear equations by implementing

$$\int\_{0}^{1} \psi\_{j}(t) R\_{i}(t) dt = 0, \; j = 1, \ldots, 2^{k} (M+1) - mn, \; i = 1, 2, \ldots, m. \tag{32}$$

Also, by substituting the initial conditions in Equation (28) into Equation (30), then we obtain

$$\begin{aligned} &u\_{1}(t\_{0}) \approx \mathbb{C}\_{1}^{\top}\boldsymbol{\Psi}(t\_{0}) = u\_{10}, \qquad \frac{d\boldsymbol{u}\_{1}}{dt}(t\_{0}) \approx \mathbb{C}\_{1}^{\top}D\boldsymbol{\psi}(t\_{0}) = u\_{11}, \qquad \dots, \quad \frac{d^{n-1}u\_{1}}{dt^{n-1}}(t\_{0}) \approx \mathbb{C}\_{1}^{\top}D^{n-1}\boldsymbol{\psi}(t\_{0}) = u\_{1(n-1)} \\ &u\_{2}(t\_{0}) \approx \mathbb{C}\_{2}^{\top}\boldsymbol{\psi}(t\_{0}) = u\_{20}, \qquad \frac{d\boldsymbol{u}\_{2}}{dt}(t\_{0}) \approx \mathbb{C}\_{2}^{\top}D\boldsymbol{\psi}(t\_{0}) = u\_{21}, \qquad \dots, \quad \frac{d^{n-1}u\_{2}}{dt^{n-1}}(t\_{0}) \approx \mathbb{C}\_{2}^{\top}D^{n-1}\boldsymbol{\psi}(t\_{0}) = u\_{2(n-1)} \\ &\vdots \quad \vdots \quad \vdots \quad \vdots \quad \vdots \quad \vdots \quad \vdots \\ &u\_{m}(t\_{0}) \approx \mathbb{C}\_{m}^{\top}\boldsymbol{\psi}(t\_{0}) = u\_{m0}, \quad \frac{d\boldsymbol{u}\_{m}}{dt}(t\_{0}) \approx \mathbb{C}\_{m}^{\top}D^{\dagger}\boldsymbol{\psi}(t\_{0}) = u\_{m1}, \qquad \dots, \quad \frac{d^{n-1}u\_{m}}{dt^{n-1}}(t\_{0}) \approx \mathbb{C}\_{m}^{\top}D^{n-1}\boldsymbol{\psi}(t\_{0}) = u\_{m(n-1)} \end{aligned} \tag{33}$$

A 2*k*(*M* + 1) set of linear equations is generated by combining Equations (32) and (33). The solution of these linear equations can be obtained for unknown coefficients of the vector *C*. Consequently, *u*1(*t*), *u*2(*t*),..., *um*(*t*), introduced in Equation (27), can be computed.

If *Ui* is a nonlinear function of *t*, *u*1, *u*2, ... , *um*, then we first compute *R*1(*t*), *R*2(*t*), ... , *Rm*(*t*) at <sup>2</sup>*k*(*<sup>M</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> *mn* points and, for a better result, use the first 2*k*(*<sup>M</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> *mn* roots of shifted Legendre *<sup>P</sup>*2*<sup>k</sup>* (*M*+1)(*t*). Then these equations, collectively with Equation (33), produce 2*k*(*<sup>M</sup>* + <sup>1</sup>) nonlinear equations. The solution of these nonlinear equations can be obtained by employing Newton's iterative method. Consequently, *u*1(*t*), *u*2(*t*),..., *um*(*t*), introduced in Equation (27), can be computed.
