Numerical Solutions for Systems of Fractional and Classical Integro-Differential Equations via Finite Integration Method Based on Shifted Chebyshev Polynomials

: In this paper, the ﬁnite integration method and the operational matrix of fractional integration are implemented based on the shifted Chebyshev polynomial. They are utilized to devise two numerical procedures for solving the systems of fractional and classical integro-differential equations. The fractional derivatives are described in the Caputo sense. The devised procedure can be successfully applied to solve the stiff system of ODEs. To demonstrate the efﬁciency, accuracy and numerical convergence order of these procedures, several experimental examples are given. As a consequence, the numerical computations illustrate that our presented procedures achieve signiﬁcant improvement in terms of accuracy with less computational cost.


Introduction
Fractional calculus is a branch of mathematical analysis that has been receiving much attention from many researchers. Due to the fact that several real-world phenomena can be described successfully by developing mathematical models using fractional derivatives and integrations. Some interesting applications of fractional calculus can be found in various fields of sciences and engineering, for examples, viscoelasticity [1], nonlinear dynamical system [2], chaotic system [3], electromagnetic wave [4], heat transfer modeling [5], etc. In addition, the fractional calculus has grown in popularity as a tool for describing the physical features of real-world situations, particularly COVID-19, SIR model and health problems, as evidenced in [6,7] and references therein. One interesting issue regarding the fractional calculus is a fractional integro-differential equation (FIDE). It consists of both integral and differential operators involving derivatives of positive fractional order. The FIDEs have demonstrated to be adequate models for several phenomena arising in damping laws, earthquake model, diffusion processes, fluid dynamics, traffic models and acoustics, see [8][9][10] and references cited therein for more details. However, the fractional order derivative of FIDEs can be reduced to a positive integer order. Then, it is called the classical integro-differential equation (CIDE) which is frequently used to describe many applications which can be seen in [11][12][13] for details of applications. Actually, many problems of both FIDE and CIDE are often constructed to be a system.
In fact, most of FIDEs and CIDEs and system involving them are usually difficult to solve analytically. Therefore, numerical techniques are required to obtain an accurate approximate solution. Several numerical methods for solving the FIDEs and CIDEs have been given, for examples, variational iteration method [8], collocation method [13], homotopy method [14], Adomian's decomposition method [15]. In 2013, an efficient numerical method had occurred which is called the finite integration method (FIM) introduced by Wen et al. [16]. It has been developed in order to solve one-dimensional partial differential equations (PDEs). The concept of FIM is to transform a given PDE into an equivalent integral equation and then numerical integrations are applied. It is known that the integration task involves multiplication by a small step size, whereas the differentiation task involves division by a small step size. As a reason, the numerical integration is very insensitive to round-off error and preserves the approximation accuracy. Therefore, the approximation of FIM can provide a stable, accurate and efficient numerical solution, see [16][17][18]. In 2015, this FIM has been extended to overcome the multi-dimensional PDEs found in [17]. After that, the FIM has been improved by hiring numerical quadratures such as Simpson's rule, Newton Cotes and Lagrange interpolation, see [18]. As a consequence, these improved FIMs give highly accurate solutions compared with the traditional FIM and finite difference method (FDM). In 2018, Boonklurb et al. [19] have modified the FIM by using Chebyshev polynomials to solve one-and two-dimensional PDEs. This modified FIM also provides much higher accuracy than the FDM and those original FIMs with small computational nodes. Recently, the modified FIM was widely utilized to apply with many applications, see [20][21][22][23]. Also, it was demonstrated that results obtained by the modified FIM achieve significant improvement in terms of accuracy more than several existing methods.
The major aim of this paper is to develop the modified FIM [19] by using the shifted Chebyshev polynomial which thereafter will be referred to as FIM-SCP and also constructs the operational matrix of fractional integration in order to devise two numerical procedures for solving numerically the systems of both FIDEs and CIDEs of the Volterra type. Actually, the technique of FIM-SCP has been proposed in [20]. It is used to find numerical solutions of direct and inverse problems for the time-dependent Volterra integro-differential equation (VIDE). Hence, in this paper, we continue our study from [20] by extending the VIDEs to be a kind of system that involves the fractional-order differential operator. The problem mainly considered in this article is called the system of FIDEs which is studied in the form presented in [24], i.e., for all i ∈ {1, 2, 3, . . . , m} and L ∈ R + with initial conditions y (v) is parameter describing the order of fractional derivative D α i in the Caputo sense [25], α i is the smallest integer greater than α i , f i (x) and p ij (x) are sufficiently continuous functions, κ ij (x, t) is continuously integrable kernel function and y i (x) is unknown function to be solved numerically. Next, when the derivative of fractional order is reduced to any α i ∈ N, we obtain the system of CIDEs. In this paper, we investigate the system of CIDEs in the following general form for all i ∈ {1, 2, 3, . . . , m} and L ∈ R + with initial conditions y specified constant and h i = max 1≤j≤m r ij when r ij is the highest order of derivative for each y j contained in the linear differential operator L ij which is defined by (33), λ ij is constant coefficient, κ ij (x, t) is continuously integrable kernel function, f i (x) is continuous function and y j (x) is unknown function to be determined. Moreover, we observe that if the kernel functions κ ij (x, t) or the constant coefficients λ ij in (2) are all zeros, (2) also becomes the system of ordinary differential equations (ODEs). An interesting problem for the system of ODEs is the stiff problem which can be seen in [26] for modeling various real-world problems. Nevertheless, the stiff system of ODEs is difficult to solve analytically and numerically. A stiff system generally happens when some components of the solutions decay much more rapidly than others. It affects their numerical solutions in terms of stability. This feature forces the used numerical method to choose an extremely small step size which consumes the computational times expensively and may give inaccurate solutions. Accordingly, examples of the stiff system of ODEs are also presented to illustrate the efficiency of the proposed procedure that can be treated these troubles. In this study, we assume that (1) and (2) have unique solutions under the given supplementary conditions. The organization of this paper is as follows. In Section 2, the developed FIM-SCP given by [20] is briefly introduced which is utilized to be the principal tool for devising the numerical procedures. In Section 3, the shifted Chebyshev expansion is employed to construct the operational matrix of fractional integration. It is together used with the FIM-SCP to devise the procedure for solving the system of FIDEs (1). This procedure is verified the efficiency with several examples. Next, Section 4 includes the procedure for solving the system of CIDEs (2) and experimental examples for testing accurate solutions obtained by the procedure. This procedure is also applied with the stiff system of ODEs via several examples in this section. Finally, the conclusion and discussion are summarized in Section 5.

The Developed FIM-SCP
In this section, we briefly introduce the technique of FIM-SCP presented in [20] which is utilized to be the principal tool for devising the numerical procedures to solve systems (1) and (2). Let us introduce the definition of shifted Chebyshev polynomials. Definition 1 ([27]). The shifted Chebyshev polynomial of degree n ≥ 0 is defined by Moreover, the analytic form of S n (x) with n > 0 given by [27] can be written as Some important properties of the shifted Chebyshev polynomial are further given in Lemma 1. They will be used to construct the first and higher orders of the shifted Chebyshev integration matrices which are the major tools of the FIM-SCP. Lemma 1 ([20]). The followings are properties of the shifted Chebyshev polynomial (3).
(i) The zeros of shifted Chebyshev polynomial S n (x) for x ∈ [0, L] are (ii) The vth order derivatives of shifted Chebyshev polynomial S n (x) at x = 0 are (iii) The single integrations of shifted Chebyshev polynomial S n (x) for n ≥ 2 are where the initial integrations of S 0 and S 1 are x and x 2 L − x, respectively.
Next, we construct the first order integration matrix by letting j and M be natural numbers, u i (x) be an approximate solution of y i (x) contained in (1) and (2) which is defined by a linear combination of the shifted Chebyshev polynomials S 0 (x), S 1 (x), S 2 (x), . . . , S M−1 (x). Then, we have where c n i is the unknown coefficient to be considered. Let x k be grid points generated by the zeros of the shifted Chebyshev polynomial S M defined by (5) in ascending order. When we substitute each x k into (8), they can be expressed in the matrix form which is denoted by u i = Sc i . Since S is invertible by Lemma 1(iv), we have c i = S −1 u i . Now, we consider the single integration of u i from 0 to x k , denoted U i (x k ), to obtain where S n is denoted to the single-layer integration of S n that can explicitly find by (7) depending on its degree n. After substituting each node x k into the above equation, it can be written in the matrix form M×M is the integral operational matrix called the first order shifted Chebyshev integration matrix (SCIM). It can be also expressed to another form Remark 1 ([20]). Based on (9), the m-layer integration of u i from 0 to x k , denoted U (m) i (x k ), can be easily obtained by the first order SCIM A multiplied by itself m times, i.e.,

Procedure for Solving System of FIDEs
In this section, we devise the numerical procedure for solving the system of FIDEs (1) and verify accurate solutions obtained from this proposed procedure by comparing with other existing methods and their analytical solutions. Before constructing the procedure, we attempt to create an operational matrix of integrals for the fractional integration and derivative. Let us provide some basic definitions and properties of fractional calculus theory from [25,28] as follows.
Definition 2. Let p, µ and x be real numbers such that x > 0 and The Riemann-Liouville fractional integral operator of order α of u ∈ C µ , µ ≥ −1 is defined by where Γ(·) is the well-known Gamma function.

Definition 3.
The Caputo fractional derivative D α of u ∈ C n −1 for n ∈ N is defined by Actually, we have known that the Riemann-Liouville fractional integral operator I α is a linear operator which has some important properties including x k+α for k ∈ N ∪ {0}. (11) Recall that for α ∈ N, the Caputo differential operator D α coincides with the usual differential operator of integer order. More properties can be found in [25,28].

Operational Matrix of Fractional Integration
In this work, the fractional derivatives of the system of FIDEs (1) are studied based on the Caputo sense stated in Definition 3. We create the the operational matrix of fractional integrals by using the shifted Chebyshev expansion which is called the shifted Chebyshev fractional matrix (SCFM). First, we express and prove the analytic formula of the Riemann-Liouville fractional integral with order α for the shifted Chebyshev polynomials as demonstrated in Theorem 1. Theorem 1. Let S n (x) be a shifted Chebyshev polynomial. Then, Proof of Theorem 1. For n = 0, it is obvious that S 0 (x) = 1. Then, using Definition 2 for α > 0, we get In addition, for n > 0, we have known that the fractional integration in the Riemann-Liouville sense is a linear operation. Thus, by employing (4) and (11), we obtain Hence, a combination of (13) and (14) leads to the desired result (12).
Next, we construct an operational matrix of the Riemann-Liouville fractional integral with order α for an approximate unknown function u i (x) by utilizing the shifted Chebyshev expansion (8) which is called the SCFM as expressed in Theorems 2 and 3.

Theorem 2.
Let u i (x) be approximated by the shifted Chebyshev expansion (8). Then, where is an M row vector in which each component can be calculated by (12), is an M column vector of each u i (x) at the zeros x k in (5) and S −1 is defined in Lemma 1(iv).

Proof of Theorem 2.
Since the Riemann-Liouville's fractional integration of order α > 0 is a linear operation and use the linear combination u i (x) in (8), we have Rewriting the above equation as a vector form and use the vector c i in Section 2, we get Hence, we obtain the desired representation (15). Theorem 3. Let u i be an M column vector of each u i (x) at the zeros x k in (5). Then, Proof of Theorem 3. Let α > 0 and M ∈ N. By employing the relation (15), we obtain Thus, we achieve the operational matrix of Riemann-Liouville fractional integral. Note that, the elements contained in J α can be computed by (12).
However, we observe that for the order α = 1, we have J α = S which is defined in Section 2. Thus, when α is a positive integer, then A α = J α S −1 . In order to reduce the computational time for the positive integer order α, we consume the SCIM A α instead of the SCFM J α S −1 . Because, when the order α has changed, the matrix J α needs to recalculate its elements again. Conversely, for the matrix A α , the SCIM A is computed only once, then it is raised to the power α. Hence, we obtain the following Remark 2.

Numerical Solution for System of FIDEs
The objective of this section is to create a numerical procedure for solving the system of FIDEs (1) with the given initial conditions. First, let u i (x) be an approximate solution of y i (x) for (1). Thus, (1) can be written in the form for all i ∈ {1, 2, 3, . . . , m} and L ∈ R + with the initial conditions First, we discretize the domain [0, L] into M grid nodes which are generated by the zeros of S M defined in (5), where x 1 < x 2 < x 3 < · · · < x M . Then, the approximate solution u i (x) is sought at these zeros. To simplify (16), we denote each integral term by for j ∈ {1, 2, 3, . . . , m}. Thus, (16) becomes Next, we attempt to eliminate the fractional derivatives from (19) by taking the α ilayer integrals from 0 to a zero x k on both sides of (19). Then, we have . . , d i α i are arbitrary constants emerged in the process of integrations. Conveniently, we can rewrite the above equation in another form by using the integral operator I instead. Then, it becomes where the integration of Caputo fractional derivative term has gotten by using (10).
After that, let us consider each term contained in (20) and reformulate it to the matrix form by varying the variable x k as the zeros x 1 , x 2 , x 3 , . . . , x M which uses the notation " −→" for representing the mapping from R to R M . Then, we have . . .
where J α i −α i and S −1 are defined in Theorem 3 and Lemma 1(iv), respectively. Next, we consider the summation term of the unknown constants in (20), that is, Then, the integration of the forcing terms in (20) is formulated by using Remark 2 to the matrix form . . .
After that, we consider the remaining terms in the summation on the right-hand-side of (20). Then, its first term is converted to the matrix form . . . where For another term in the summation on the right-hand-side of (20), I α i Q ij (x k ), before transforming it to the matrix form, we consider Q ij (x k ) by using (9) and (18) to obtain where a kl is an element at the kth row and the lth column of the SCIM A = SS −1 defined in Section 2. Thus, we have the column vector of Q ij (x k ) as is the Hadamard product defined in [29].
From the above relation and Remark 2, we then have Consequently, by substituting expressions (21)- (25) into (20), we obtain or it can be simplified as where H ij := P ij + (A K ij ). Finally, we vary all indices i ∈ {1, 2, 3, . . . , m} in (26). Then, we have the following system which the system can be rearranged to the block-matrix of the form where I m is an m × m identity matrix, S −1 is defined in Lemma 1(iv) and ⊗ is the Kronecker product defined in [29]. Other parameters contained in (27) are defined by the following block matrices: where "blkdiag · " is a block diagonal matrix in which the off-diagonal elements are the zero matrices. However, we can see that (27) has unknown vectors apart from u, i.e., d which is emerged from the process of integration for a total of ∑ m i=1 α i elements. Therefore, we require ∑ m i=1 α i equations more which are constructed by using the given initial conditions (17). At specified index i, we use (8) to transform these conditions into the matrix form, we have u M−1 (0) is the M row vector in which its elements can be found by using (6), . . , u i (x M )] and S −1 is defined in Lemma 1(iv). Thus, when all derivative orders v ∈ {0, 1, 2, . . . , α i − 1} are varied in the above equation, we obtain (1) . . .
Then, we substitute all indices i ∈ {1, 2, 3, . . . , m} in (28) and write them into the block-matrix form where I m is an m × m identity matrix, S := blkdiag S 1 , S 2 , S 3 , . . . , S m and b = b 1 , b 2 , b 3 , . . . , b m . Since (28) has α i equations, it now implies that (29) has ∑ m i=1 α i equations. Therefore, we achieve more equations as required which resulted in the numbers of unknown variables and the numbers of equations in (27) and (29) to be precisely equivalent. Finally, we can combine both (27) and (29) in order to construct the linear system in the block-matrix form as the following where 0 is the ∑ m i=1 α i × ∑ m i=1 α i zero matrix and each parameter contained in (30) is defined as mentioned above. Note that the amount of calculation for the linear system (30) totally consists of mM + ∑ m i=1 α i equations which can be certainly solved by using the backslash or inverse command in MatLab solver. Now, we can seek the approximate solution u by solving (30). The obtained solution u consists of each result u i for i ∈ {1, 2, 3, . . . , m} which is at the position followed the zeros of shifted Chebyshev polynomial S M . Nevertheless, if we would like to know the solution u i (x) at a different position within the domain [0, L], it can be calculated by where S(x) := S 0 (x), S 1 (x), S 2 (x), . . . , S M−1 (x) and u i obtains from solving (30).

Experimental Examples for System of FIDEs
In order to illustrate the effectiveness of the proposed numerical procedure in the preceding section, we now present some experimental examples for solving the system of FIDEs (1). In calculation, we implement the proposed method to solve four examples and show its accuracy and efficiency which is measured by the absolute error

Example 1 ([24]
). Consider the following linear system of FIDEs over x ∈ [0, 1] for α 1 , α 2 ∈ (0, 1] with the initial conditions u 1 (0) = 1 and u 2 (0) = −1. The exact solutions for In this Example 1, the fractional orders of derivative are considered on (0, 1]. By employing our numerical procedure, we have the approximate solutions u 1 (x) and u 2 (x) for any x ∈ [0, 1]. In Table 1, the absolute errors between the exact and approximate solutions at α 1 = α 2 = 1 are demonstrated. It also shows a comparison of the absolute errors with M = 16 between our suggested procedure and the technique based on operational matrices of triangular functions (OMTF) given by [24]. We can see that our method provides much higher accuracy than the OMTF. The consuming time for M = 16 via MatLab program is about 0.0805 s. Additionally, Figure 1a,b depict the comparisons between analytical and numerical solutions u 1 (x) and u 2 (x) with M = 30 and α 1 = α 2 = 1. We can see that our obtained solutions quite match exactly. Table 1. Absolute errors of u 1 (x) and u 2 (x) when α 1 = α 2 = 1 for Example 1.
x In Table 2, absolute errors of the approximate solutions u 1 (x) and u 2 (x) are shown at different orders α 1 = α 2 := α ∈ {0.99, 0.999, 0.9999} and M = 16. This investigation shows that as the fractional order α increases from 0.99 to 0.9999, the respective accuracy is increasing and attained its maximum accuracy at α = 1. Finally, the error analysis again verifies that at several values of α ∈ {0.91, 0.93, 0.95, 0.97, 0.99, 1}, the obtained numerical solutions converge to the integer order solutions which their behaviors are also displayed in Figure 1c,d. They indeed attain to the blue solid line (α = 1).

Example 2 ([30]
). Consider the following linear system of FIDEs over x ∈ [0, 1] for α 1 , α 2 ∈ (0, 2] with the initial conditions u 1 (0) = 1, u 2 (0) = 0, u 1 (0) = 1 and u 2 (0) = 2. The analytical solutions when α 1 = α 2 = 2 are u * 1 (x) = x + cos x and u * 2 (x) = x + sin x.  This experimental example is considering the problem under the fractional orders of derivatives within the interval (0, 2]. By using our proposed numerical procedure, we have the numerical solutions u 1 (x) and u 2 (x) for each x ∈ [0, 1] which are compared with solutions given by the OMTF [24] measured by the absolute error as shown in Table 3. We can see that our FIM-SCP produces higher accuracy than the OMTF with the same number of discretization nodes M = 16 when the fractional orders α 1 = α 2 = 2. Moreover, we plot the obtained approximate solutions u 1 (x) and u 2 (x) compared to their analytical solutions as depicted in Figure 2a,b, respectively. These figures show that these solutions perfectly match with the exact values. The computational time for M = 16 of this process is around 0.0984 s. Table 3. Absolute errors of u 1 (x) and u 2 (x) when α 1 = α 2 = 2 for Example 2.
x Furthermore, Table 4 demonstrates the absolute errors of approximate solutions u 1 and u 2 at different values of the fractional orders α 1 = α 2 := α ∈ {1.99, 1.999, 1.9999} with M = 16. We can see that their errors are decreasing when α → 2. Finally, Figure 2c 1.93, 1.95, 1.97, 1.99, 2}. We can see that they also tend to the blue solid line, that is, the integer order solution at α = 2. Table 4. Absolute errors of u 1 (x) and u 2 (x) at different orders α for Example 2.

Example 3 ([30]
). Consider the following linear system of FIDEs over x ∈ [0, 1] This system consists of three equations of FIDEs and the fractional orders are in the interval (0, 1]. Based on the presented numerical procedure, by solving the corresponding block-matrix Equation (30) with M = 16 and α 1 = α 2 = α 3 = 1, we obtain the numerical solutions u 1 (x), u 2 (x) and u 3 (x). When we find accuracy of these obtained solutions via the absolute error and compare with the OMTF [24], we can see that our method gives a much higher accuracy as shown in Table 5  From Examples 1 (α 1 = α 2 = 1), 2 (α 1 = α 2 = 2) and 3 (α 1 = α 2 = α 3 = 1), we can see that our proposed procedure for solving the system of FIDEs provides an excellent accuracy in terms of the absolute error when the fractional orders tend to integer orders which can be seen in the previous tables. Moreover, we can say that a sequence (u M ) converges to the exact solution u * with order p if there exists a constant C such that u * − u M < CM −p or u * − u M = O(M −p ) using the big-O notation [31]. In practice, for approximating an order of convergence p, we take the natural logarithmic function on both sides of the above expression. Thus, we obtain ln u * − u M ≈ ln C − p ln M. However, if we take two distinct discretizations M old and M new into the equation, we can solve them to find the estimation of convergence order by p ≈ , where u M old and u M new denote numerical solutions obtained by using the consecutive discretization nodes M old and M new in ascending order, respectively. Hence, the orders of convergence for Examples 1-3 based on Euclidean norm are considered numerically and presented in Table 6. It is obvious that the obtained convergence orders rapidly increase when the number of nodes M ever increases.  Table 5. Absolute errors of u 1 (x), u 2 (x) and u 3 (x) at α 1 = α 2 = α 3 = 1 for Example 3. x

Procedure for Solving System of CIDEs
In this section, we extend the concept of solving system of FIDEs in Section 3 by studying in the common case, i.e., the order of the fractional derivative is focused on the positive integer which is called CIDE. Therefore, we construct the numerical procedure for solving the generalized system of CIDEs (2). An accuracy of the obtained solutions is also verified in this section.

Numerical Solution for System of CIDEs
Let us first introduce the system of CIDEs (2) with the given initial conditions. Let u j (x) be an approximate solution of y j (x) contained in (2). Then, it becomes for all i ∈ {1, 2, 3, . . . , m} and L ∈ R + with the given initial conditions u where h i = max 1≤j≤m r ij when r ij is the highest order of derivative for each u j contained in the linear differential operator L ij which is defined by where p r ij (x) for each r ∈ {0, 1, 2, . . . , r ij } is continuously differentiable function up to the highest order of derivative contained in (32). We apply the idea of FIM-SCP described in Section 2 to deal with the integration term in (32). Then, the numerical procedure for solving the system of CIDEs is constructed. First, by using (18), we can rewrite (32) into Next, we discretize the domain [0, L] into M nodes which are generated by the zeros of S M defined in (5), i.e., x 1 < x 2 < x 3 < · · · < x M . Using the FIM-SCP, we have to eliminate all derivatives from (34). Actually, we know that the highest order of derivatives in each ith equation of (34) that is h i . In order to remove its derivatives, the h i -layer integrals from 0 to x k are taken on both sides of (34). Thus, we obtain where d i1 , d i2 , d i3 , . . . , d ih i are arbitrary constants emerged in the process of integrations and Z ij (x k ) is an integration terms for each u j . By employing (33), it can be defined by where p r ij is the coefficient function corresponding to the rth order derivative of u j for r ∈ {0, 1, 2, . . . , r ij }. Next, we reformulate Z ij (x k ) using the integration by parts for each term to remove all derivatives from u j . Anywise, Z ij (x k ) can be separated into two cases which are considering as r ij = h i and r ij < h i . Thus, for the first case r ij = h i , it becomes and for the second case r ij < h i , we have where (p r ij ) (n) is the nth order derivative of p r ij for r, n ∈ {0, 1, 2, . . . , r ij } in both (36) and (37). Hence, (35) can be written in another form where I h i is the h i -layer repeated integral operator from 0 to the zero x k . Subsequently, we apply the idea of our proposed FIM-SCP described in Section 2 to transform (38) into the matrix form by varying x 1 , x 2 , x 3 , . . . , x M as the zeros of shifted Chebyshev polynomial S M defined in (5).
Then, for the first case r ij = h i , we can express (36) in the matrix form and for the second case r ij < h i , (37) can be written in the matrix form ] . However, we can explicitly see that when substituting r ij = h i in (40), it indeed becomes (39). Thus, we can only use (40) which is enough to represent Z ij for both cases, i.e., r ij ≤ h i . Moreover, (40) can be simplified to Next, we transform the remaining terms of (38) into the matrix form by utilizing the same processes with (22), (23) and (25) that change α i into h i instead, respectively, i.e., where the parameters are defined in the similar idea as presented in Section 3.2. Hence, by substituting the expressions (41) and (42) into (38), we obtain . Finally, when it is varied as all indices i ∈ {1, 2, 3, . . . , m}, we have where the parameters contained in (43), except H, are defined as same as (27) in which change α i into h i instead and For the given initial conditions, we can perform as same processes as to obtain (28) and (29) in the previous section. We finally have S I m ⊗ S −1 u = b, where the parameters are defined as in Section 3.2 by instead changing α i to h i . Note that, it has exactly ∑ m i=0 h i equations. In practice, we can solve the systems (43) and above condition, simultaneously, by constructing them to the linear system in a block matrix form where 0 is the ∑ m i=1 h i × ∑ m i=1 h i zero matrix and other parameters are defined as mentioned above. Note that the linear system (44) has a total of mM + ∑ m i=1 h i equations. Hence, by solving (44), we obtain the numerical solution u i (x) for the problem (2).

Experimental Examples for System of CIDEs
In this section, we apply the proposed numerical procedure to find approximate solutions of the systems of CIDEs (2). We numerically demonstrate three examples that are computed via MatLab program to verify the accuracy of our algorithm.

Example 4 ([32]
). Consider the system of linear first order CIDEs over x ∈ [0, 1] subject to the initial conditions u 1 (0) = 1 and u 2 (0) = −1. Whose analytical solutions of this example are given by u * 1 (x) = x + e x and u * This experiment is the system of first order CIDEs which is considered under constant coefficients, constant kernel functions and polynomial forcing terms. By using our numerical procedure, we obtain the approximate solutions u 1 (x) and u 2 (x) for each x ∈ [0, 1]. A comparison of the absolute error between u 1 (x) and u 2 (x) obtained from our proposed method, Genocchi polynomials method (GPM) [32] and biorthogonal system in approximation (BSA) [33], with their exact solutions by using M = 8 as shown in Table 7. The run-time is about 0.0437 s.
This example is a system of linear second order CIDEs with variable coefficients, polynomial forcing terms and kernel functions are in term of functions depending on variables x and t. By hiring our procedure with M = 8, we obtain the numerical solutions u 1 (x) and u 2 (x) for each x ∈ [0, 1]. We compare the absolute errors between our algorithm and STWS [34] as displayed in Table 8. The run-time is 0.0880 s. From Examples 4 and 5, we can see that our proposed method for solving the system of CIDEs provides higher accuracy than other existing methods in terms of the absolute errors at the same number of nodes and under the same conditions which can be seen in Tables 7 and 8. Their average run-times are also consumed very inexpensive. In addition, we demonstrate the Euclidean error norm and the convergence orders p or O(M −p ) of the obtained solutions from Examples 4 and 5 in Table 9 for varying nodes M. We can see that our proposed procedure acquires a significant improvement in term of accuracy with less computational nodes M. This table also shows that the convergence orders rapidly increase when the number of nodes M ever increases. Actually, we can also apply the devised numerical procedure for solving system of CIDEs to overcome the stiff system of ODEs by vanishing the kernel functions κ ij or the parameters λ ij . Thus, we further illustrate two examples of the stiff systems. Example 6. Consider the following stiff system of linear first order ODEs over x ∈ [0, 1]. u 1 (x) = 98u 1 (x) + 198u 2 (x), u 2 (x) = −99u 1 (x) − 199u 2 (x), subject to the initial conditions u 1 (0) = u 2 (0) = 1. The analytical solutions are u * 1 (x) = 4e −x − 3e −100x and u * 2 (x) = −2e −x + 3e −100x .
By employing the proposed algorithm, we obtain the approximate solutions u 1 (x) and u 2 (x) at different x close to zero as shown in Table 10. When they are compared with the Runge-Kutta fourth-order method (RK4) at M = 16, we have that our obtained results give higher accuracy. We also capture the solutions with M = 30 near zero as depicted in Figure 4a,b. Moreover, we found that the RK4 method is unavailable for solving this problem whole domain [0, 1] when small discretization. In contrast, our procedure can treat this trouble as illustrated in Figure 4c,d with M = 50. It is shown that these found solutions perfectly match with the exact.

Conclusions
In this article, the FIM-SCP is implemented to devise two numerical procedures for solving the systems of linear FIDEs and CIDEs. For the system of FIDEs (1), the fractional derivative is considered in the Caputo sense which is manipulated by the novel operational matrix of fractional integration (SCFM) as shown in Theorem 3. Then, the first procedure is devised by combining FIM-SCP and SCFM to overcome the systems of FIDEs (1) as expressed by the linear system (30). Next, the second procedure is created to deal with the generalized system of CIDEs (2) based on the FIM-SCP as expressed by the linear system (44). Note that both proposed numerical procedures are based on the FIM-SCP which completely eliminate the dilemma in the well-known round-off and discretization errors. These procedures are also in form of the linear systems. Thus, it is very convenient to find approximate solutions, we just substitute parameters of a given problem and supplementary conditions into those linear systems and solve them by MatLab solver. Moreover, we can see from several Examples 1-5 that the devised algorithms provide much higher accuracy than other methods, consume low computational cost in terms of CPU time(s) and acquire a significant improvement in terms of absolute error with less computational nodes M. The obtained convergence orders for each example also give highly order. In addition, the second approach can solve the stiff system of linear ODEs as shown in Example 6. The obtained solutions are in good agreement with the exact solution. An interesting direction for our future work is to extend our techniques to solve the multi-dimensional FIDEs and CIDEs.