**1. Introduction**

In recent years, fractional calculus has becoming an efficient and successful tool for the analysis of several physical-mathematical problems. The main reason for the increasing number of papers dealing with fractional problems is also explained by the intrinsic and natural possibility of the fractional calculus to take into account also some memory effects, which is quite impossible by using the ordinary differential operators [1]. In this work, we consider the nonlinear multi-order fractional differential equations (MOFDEs) of the form

$$\mathcal{D}^{\gamma}\_{\star} \mathbf{x}(t) = F\left(t, \mathbf{x}(t), \mathcal{D}^{\beta\_1}\_{\star} \mathbf{x}(t), \dots, \mathcal{D}^{\beta\_\ell}\_{\star} \mathbf{x}(t)\right), \quad 0 \prec t \lessapprox L,\tag{1}$$

subjected by the following boundary or supplementary conditions

$$H\_j(\mathbf{x}(\eta\_j), \mathbf{x}^{(1)}(\eta\_j), \dots, \mathbf{x}^{(m-1)}(\eta\_j)) = d\_j \quad j = 0, 1, \dots, m - 1,\tag{2}$$

where *Hj* are linear functions, *dj* P R, and *ηj* are some given points in r0, *<sup>L</sup>*s. In (1), <sup>D</sup>*γ*‹ denotes the Caputo fractional derivative operator such that *m* ´1 ă *γ* ď *m*, *m* P N, and 0 ă *β*1 ă *β*2 ă ... ă *β* ă *γ* are real constants. The function *F* can be either linear or nonlinear function of its arguments. In [2], some preliminary results both on the existence and uniqueness of the solution of MOFDEs (1) are obtained.

It is well-known that usually the exact solution of fractional differential equations cannot be obtained analytically. Therefore, many authors have recently developed some suitable numerical methods for such equations. Among the many approximation algorithms for (1) and (2), we mention the systems-based decomposition approach [3], the Adomian decomposition method [4], the spectral methods [5–8], the B-spline approach [9], and the generalized triangular function [10].

It is known that the traditional orthogonal polynomials such as Jacobi, Hermit, and Laguerre are solution of a second order differential equation. In addition, the derivatives of these polynomials constitute an orthogonal system. Moreover, there exist another set of polynomials with the two aforementioned properties. They satisfy the following second order differential equation

$$\left(\mathbf{x}^2 \, y''(\mathbf{x}) + 2(\mathbf{x} + 1)\, y'(\mathbf{x}) - n(n+1)\, y(\mathbf{x}) = 0,\tag{3}$$

where *n* is a positive integer. Krall and Frink [11] called these the Bessel polynomials thank to their close relation with the Bessel functions of half-integral order. In fact, they have shown that these polynomials occur naturally as the solutions of the classical wave equation in spherical coordinates. These polynomials also appear in the study of various mathematical topics including transcendental number theory [12,13] and student t-distributions [14]. These polynomials seem to have been considered first by Bochner [15] as pointed in Grosswald [16]. However, Krall and Frink considered them in more general setting so that to include also the polynomial solutions of the differential Equation (3). The properties of Bessel polynomials such as recurrence relations, generating functions, and orthogonality were investigated in [11]. The algebraic properties of these polynomials were considered by Grosswald [16]. Some more information about the theory and applications of Bessel polynomials can be found e.g., in [17].

In this research, we establish a new approximation algorithm based upon the Bessel polynomials to obtain a solution of a fractional model (1). In fact, one of our motivation comes from a recent paper [18], which proved the total positiveness of any collocation matrix of theses polynomials evaluated at positive (collocation) points. To the best of our knowledge, this is the first attempt to study these polynomials for approximating MOFDEs. In summary, the main idea behind the presented approximation algorithm based on using the Bessel polynomials with together the collocation points is that it transforms the differential operators in (1) to an equivalent algebraic form, thus greatly reducing the numerical efforts. It should be mentioned that our Bessel polynomials are different from those Bessel functions known as Bessel functions of the first kind that previously considered in the literature, see [19] for a recent review.

The content of the paper is structured as follows. In Section 2 some relevant properties of the Caputo fractional derivative and the generalized Taylor's formula for the Caputo derivative are presented. Section 3 is dedicated to the definitions of Bessel polynomials and their generalized fractional-order counterpart. Moreover, the results about the convergence and error bound of these polynomials are established. In Section 4, where a collocation method also shown to solve the MOFDEs. By using these Bessel basis functions along with collocation points, the proposed method converts the MOFDEs into a nonlinear matrix equation. Hence, the residual error function is introduced to assess the accuracy of Bessel-collocation scheme when the exact solutions are not available. In Section 5, some examples with various parameters together with error evaluation are given to show the utility and applicability of the method. The obtained results are interpreted through tables and figures. Finally, in Section 6, the report ends with a summary and conclusion.

#### **2. Some Preliminaries**

To continue, some definitions and theorems from fractional calculus theory are presented. **Definition 1.** *Let f*p*t*q *be a m-times continuously differentiable function. The fractional derivative* <sup>D</sup>*q*‹ *of f*p*t*q *of order q* ą 0 *in the Caputo's sense is defined as*

$$\mathcal{D}\_{\star}^{q}f(t) = \begin{cases} \mathcal{Z}^{m-q} f^{(m)}(t) & \text{if } m-1 < q < m, \\ f^{(m)}(t), & \text{if } q = m, \quad m \in \mathbb{N}, \end{cases} \tag{4}$$

*where*

$$\mathcal{X}^{\emptyset}f(t) = \frac{1}{\Gamma(q)} \int\_0^t \frac{f(s)}{(t-s)^{1-q}} \, ds, \quad t > 0.$$

The properties of the operator <sup>D</sup>*q*‹ can be found in [1]. Besides the linearity, the following properties will be also used

$$\mathcal{D}\_{\star}^{\mathcal{J}}(\mathbb{C}) = 0 \quad \text{( $\mathbb{C}$  is a constant)}, \tag{5}$$

$$\mathcal{D}\_{\star}^{\mathcal{J}} t^{\beta} = \begin{cases} \frac{\Gamma(\beta + 1)}{\Gamma(\beta + 1 - q)} t^{\beta - q}, & \text{for } \beta \in \mathbb{N}\_{0} \text{ and } \beta \gg [q], \text{ or } \beta \notin \mathbb{N}\_{0} \text{ and } \beta > [q], \\ 0, & \text{for } \beta \in \mathbb{N}\_{0} \text{ and } \beta < [q]. \end{cases} \tag{6}$$

Now, we define a generalization of Taylor's formula which involves Caputo fractional derivatives (see a proof in [20]).

**Theorem 1** (Generalized Taylor's formula)**.** *Assuming that* <sup>D</sup>*k<sup>α</sup>*‹ *g*p*x*q P Cp0, *<sup>L</sup>*s*, where k* " 0, 1, ... , *N,* 0 ă *α* ď 1*, and L* ą 0*. Then, there exists a* 0 ă *θ* ď *x such that*

$$\mathcal{g}(\mathbf{x}) = \sum\_{j=0}^{N-1} \frac{\mathbf{x}^{j\mathbf{a}}}{\Gamma(j\alpha + 1)} \mathcal{D}\_{\star}^{j\mathbf{a}} \mathcal{g}(0^{+}) + \frac{\mathbf{x}^{N\mathbf{a}}}{\Gamma(N\alpha + 1)} \mathcal{D}\_{\star}^{N\mathbf{a}} \mathcal{g}(\theta), \quad \forall \mathbf{x} \in [0, L].$$

*Also, we have*

$$\left| \lg(\mathfrak{x}) - \sum\_{j=0}^{N-1} \frac{\mathfrak{x}^{j\alpha}}{\Gamma(j\alpha + 1)} \mathcal{D}\_{\star}^{j\alpha} \mathfrak{g}(0^{+}) \right| \leqslant \frac{\mathfrak{x}^{N\alpha}}{\Gamma(N\alpha + 1)} M^{\alpha} \mathcal{A}$$

*where* |D*Nα* ‹ *g*p*θ*q| ď *Mα and* D*Nα* ‹ " D*<sup>α</sup>*‹ ¨ D*<sup>α</sup>*‹ ¨¨¨ D*<sup>α</sup>*‹ *(N-times).*

Finally, we define the concept of the weighted norm used in the proof of Theorem 2:

**Definition 2.** *Let assume that g* P Cp0, *L*s *and w*p*t*q *is a weight function. Then*

$$\|\mathfrak{g}(t)\|\_{w} = \left(\int\_{0}^{L} |\mathfrak{g}(t)|^{2} \, w(t) dt\right)^{\frac{1}{2}}.$$

#### **3. Fractional-Order Bessel Functions**

In this section, definitions of Bessel polynomials as well as their generalized fractional-order version are introduced. Hence, some properties and convergence results for them are established.

#### *3.1. Bessel Polynomials*

The Bessel polynomial B*n*p*x*q of degree *n* and with constant term equal to 1 satisfies the following differential equation

$$\left(\mathbf{x}^2 \, y''(\mathbf{x}) + 2(\mathbf{x} + 1)\, y'(\mathbf{x}) - n(n+1)\, y(\mathbf{x}) = 0, \quad n = 0, 1, \dots \right)$$

Starting with B0p*x*q " 1 and B1p*x*q " 1 ` *x*, the three-terms recurrence relation for the Bessel Polynomial is

$$\mathbb{B}\_{n+1}(\mathbf{x}) = (2n+1)\mathbf{x}\,\mathbb{B}\_n(\mathbf{x}) + \mathbb{B}\_{n-1}(\mathbf{x}), \quad n = 1, 2, \dots \tag{7}$$

Beside B0p*x*q and B1p*x*q, the next four of these polynomials are listed as follows

$$\begin{aligned} \mathbb{B}\_2(\mathbf{x}) &= 1 + 3\mathbf{x} + 3\mathbf{x}^2, \\ \mathbb{B}\_3(\mathbf{x}) &= 1 + 6\mathbf{x} + 15\mathbf{x}^2 + 15\mathbf{x}^3, \\ \mathbb{B}\_4(\mathbf{x}) &= 1 + 10\mathbf{x} + 45\mathbf{x}^2 + 105\mathbf{x}^3 + 105\mathbf{x}^4, \\ \mathbb{B}\_5(\mathbf{x}) &= 1 + 15\mathbf{x} + 105\mathbf{x}^2 + 420\mathbf{x}^3 + 945\mathbf{x}^4 + 945\mathbf{x}^5. \end{aligned}$$

The coefficients of these polynomials are positive with B*n*p0q " 1 and B<sup>1</sup> *n*p0q " *n*p*n* ` 1q{2. The explicit expression for the Bessel polynomials as the unique solution of the given differential equation is defined by

$$\mathbb{B}\_n(\mathbf{x}) = \sum\_{k=0}^n \frac{1}{k!} \frac{(n+k)!}{(n-k)!} \left(\frac{\mathbf{x}}{2}\right)^k, \quad n = 0, 1, \ldots \tag{8}$$

These polynomials form an orthogonal system with respect to the weight function *w*p*x*q " expp´2{*x*q on the unite circle *C*, i.e.,

$$\frac{1}{2\pi i} \int\_{\mathcal{C}} \mathbb{B}\_{\text{fl}}(\mathbf{x}) \, \mathbb{B}\_{\text{fl}}(\mathbf{x}) \, w(\mathbf{x}) d\mathbf{x} = \frac{2(-1)^{n+1} \delta\_{nm}}{2n+1},\tag{9}$$

where *δnm* is the Kronecker delta function. Please note that the path of integration is not unique, and it can be replaced by an arbitrary curve surrounding *x* " 0. The same conclusion is true for the weight function *<sup>w</sup>*p*x*q. This implies that an arbitrary analytic function may be added to *w*p*x*q and *w*p*x*q may be multiplied by a nonzero constant. By means of the orthogonality relation (9), one may expand a function *g*p*x*q in terms of Bessel functions

$$g(\mathbf{x}) \approx \sum\_{n=0}^{\infty} a\_n \mathbb{B}\_n(\mathbf{x}),$$

where the coefficients *an* are

$$a\_{\rm II} = (-1)^{n+1} \left( n + \frac{1}{2} \right) \int\_{\mathcal{C}} \mathbb{B}\_{\rm II}(\mathbf{x}) \, g(\mathbf{x}) \, w(\mathbf{x}) d\mathbf{x}.$$

#### *3.2. Fractional Bessel Polynomials*

The fractional-order Bessel functions can be defined by introducing the change of variable *x* " *t <sup>α</sup>*{*<sup>L</sup>*, *L*, *α* ą 0 in (8). Let these polynomials will be denoted by B*α n*p*t*q " B*n*p*x*q. By generalizing these polynomials on the interval r0, *L*s we obtain

$$\mathbb{E}\mathbb{B}\_n^{\mathbf{a}}(t) = \sum\_{k=0}^n \frac{\eta^k (n+k)!}{k! \, (n-k)!} t^{\mathbf{a}\mathbf{a}}, \quad 0 \preccurlyeq t \preccurlyeq L \preccurlyeq \infty,\tag{10}$$

where *η* " 1 2*L* . It is not difficult to show that the set of fractional polynomial functions tB*α* 0,B*<sup>α</sup>* 1, ...u is orthogonal on r0, *L*s with respect to the weight function *w<sup>α</sup> L*p*t*q " *t <sup>α</sup>*´1 expp´2*<sup>L</sup>*{*<sup>t</sup> <sup>α</sup>*q. The fractional-order polynomials are useful in particular when the solutions of the underlying MOFDEs have fractional behavior.

#### *3.3. Function Approximation and Convergence*

Our goal is to obtain an approximate solution for the model problem (1) represented by the truncated Bessel series

$$\propto\_{N,a}(t) = \sum\_{n=0}^{N} a\_n \mathbb{B}\_n^a(t), \quad 0 \preccurlyeq t \preccurlyeq L,\tag{11}$$

where the unknown coefficients *an*, *n* " 0, 1, ... , *N* must be sought. For this purpose, we express B*<sup>α</sup>n*p*t*q in the matrix representation as

$$\mathbf{B}\_{\mathbf{a}}(t) = \mathbf{T}\_{\mathbf{a}}(t)\,\mathbf{D}^{t},\tag{12}$$

.

where

$$\mathbf{T}\_a(t) = \begin{bmatrix} 1 & t^a & t^{2a} & \dots & t^{Na} \end{bmatrix}, \quad \mathbf{B}\_a(t) = \begin{bmatrix} \mathbb{B}\_0^a(t) & \mathbb{B}\_1^a(t) & \dots & \mathbb{B}\_N^a(t) \end{bmatrix},$$

and the lower triangular matrix *D* of size p*N* ` 1qˆp*N* ` 1q takes the form


By expressing the relation (11) in a matrix form and exploiting (12), the approximate solution *xN*,*<sup>α</sup>*p*t*q in the matrix form can be rewritten as

$$\propto\_{\mathcal{N},a}(t) = \mathbf{B}\_a(t)\,\mathbf{A} = \mathbf{T}\_a(t)\,\mathbf{D}^t\,\mathbf{A},\tag{13}$$

where the vector of unknown is *A* " r*<sup>a</sup>*0 *a*1 ... *aN*s*<sup>t</sup>*. Our further aim is to establish the convergence results of the fractional Bessel polynomials. Roughly speaking, the next theorem shows that the approximate solution *xN*,*<sup>α</sup>*p*t*q converges to the solution *x*p*t*q of differential Equation (1) as *N* Ñ 8, see e.g., [21] for a similar proof.

**Theorem 2.** *Let assume that* <sup>D</sup>*k<sup>α</sup>*‹ *g*p*t*q P Cp0, *L*s *for k* " 0, 1, . . . , *N and let*

$$
\mathbb{S}\_N^\kappa = \operatorname{Span} \langle \mathbb{B}\_0^\kappa(t), \mathbb{B}\_1^\kappa(t), \dots, \mathbb{B}\_{N-1}^\kappa(t) \rangle.
$$

*Suppose that gN*,*<sup>α</sup>*p*t*q " *<sup>B</sup>α*p*t*q *A is the best approximation out of* S*<sup>α</sup>N to g, then the following error bound holds:*

$$\|\|\mathcal{g}(t) - \mathcal{g}\_{N,a}(t)\|\|\_{w\_L^{a}} \lesssim\_{\varepsilon} \frac{L^{Na} M^{\alpha}}{\exp(\frac{1}{L^{a-1}})\,\Gamma(N\alpha + 1)} \left(\frac{L^{a}}{2Na + a}\right)^{1/2}.$$

*where Mα* ě |D*Nα* ‹ *g*p*t*q|*, t* P p0, *<sup>L</sup>*s*.* **Proof.** According to Theorem 1, the generalized Taylor's formula for *g*p*t*q can be represented as *G* " ř*<sup>N</sup>*´<sup>1</sup> *j*"0 *tjα* <sup>Γ</sup>p*jα*`<sup>1</sup>qD*j<sup>α</sup>*‹ *<sup>g</sup>*p0`q, and satisfies

$$|\mathfrak{g} - G| \ll \frac{t^{N\alpha}}{\Gamma(N\alpha + 1)} M^{\alpha}.$$

Using the fact that *<sup>B</sup>α*p*t*q *A* is the best approximation to *g* from S*<sup>α</sup>N*and *G* P S*<sup>α</sup>N*, we conclude that

$$\left\|\boldsymbol{g}(t) - \boldsymbol{g}\_{N,a}(t)\right\|\_{\boldsymbol{w}^{\boldsymbol{\theta}}\_{L}}^{2} \leqslant \left\|\boldsymbol{g} - \boldsymbol{G}\right\|\_{\boldsymbol{w}^{\boldsymbol{\theta}}\_{L}}^{2} \leqslant \left[\frac{M^{a}}{\Gamma(Na+1)}\right]^{2} \int\_{0}^{L} \exp(-\frac{2L}{t^{a}}) t^{2Na} t^{a-1} dt.\tag{14}$$

Employing the inequality ´2*L<sup>α</sup> tα* ď ´2, which holds for all *t* P p0, *<sup>L</sup>*s, one immediately find that expp´<sup>2</sup>*Ltα* q ď expp ´2 *L<sup>α</sup>*´<sup>1</sup> q. Thus, by inserting this inequality into (14) and then integrating we conclude that

$$\|\mathcal{g}(t) - \mathcal{g}\_{N,\mathfrak{a}}(t)\|\_{\mathcal{w}\_L^{\mathfrak{a}}}^2 \lesssim \left[\frac{M^{\mathfrak{a}}}{\Gamma(N\mathfrak{a} + 1)}\right]^2 \frac{\exp(\frac{-2}{L^{\mathfrak{a}} + 1}) L^{(2N + 1)\mathfrak{a}}}{(2N + 1)\mathfrak{a}}.$$

The proof is complete by taking the square roots of both sides.

Therefore, for obtaining an approximate solution of the form (11) for the solution of (1) the following collocation points are used on 0 ă *t* ď *L*,

$$\mathbf{t}\_{i} = \frac{L}{N}\dot{\mathbf{t}}\_{\prime} \quad \dot{\mathbf{t}} = \mathbf{0}, \mathbf{1}, \dots, N. \tag{15}$$

.

#### **4. The Collocation Scheme**

To proceed, we approximate the solution *x*p*t*q of MOFDEs (1) in terms of p*N* ` <sup>1</sup>q-terms Bessel polynomials series denoted by *xN*,*<sup>α</sup>*p*t*q on the interval r0, *<sup>L</sup>*s. In the matrix representation, we consider

$$\mathbf{x}(t) \cong \mathbf{x}\_{N,a}(t) = \mathbf{T}\_a(t) \mathbf{D}^t \mathbf{A}. \tag{16}$$

By placing the collocation points (15) into (16), we ge<sup>t</sup> to a system of matrix equations as

$$\propto\_{N,a}(t\_i) = \mathbf{T}\_a(t\_i) \, \mathbf{D}^t \, \mathbf{A}, \quad i = 0, 1, \ldots, N.$$

Hence, we write the preceding equations compactly as

$$\mathbf{X} = \mathbf{T} \, \mathbf{D}^t \, \mathbf{A},\tag{17}$$

where

$$\mathbf{T} = \begin{bmatrix} \mathbf{T}\_{\alpha}(t\_0) \\ \mathbf{T}\_{\alpha}(t\_1) \\ \vdots \\ \mathbf{T}\_{\alpha}(t\_N) \end{bmatrix}, \quad \mathbf{X} = \begin{bmatrix} \mathbf{x}\_{N,\alpha}(t\_0) \\ \mathbf{x}\_{N,\alpha}(t\_1) \\ \vdots \\ \mathbf{x}\_{N,\alpha}(t\_N) \end{bmatrix}.$$

To handle the fractional derivative of order *γ* in (1), we differentiate both sides of (16),

$$\mathcal{D}\_{\star}^{\gamma} \mathbf{x}\_{\mathcal{N},a}(t) = \mathcal{D}\_{\star}^{\gamma} \mathbf{T}\_a(t) \mathbf{D}^{\dagger} \mathbf{A}. \tag{18}$$

By means of the property (5) and (6), the calculation of <sup>D</sup>*γ*‹ *<sup>T</sup>α*p*t*q can be easily obtained as follows

$$\mathbf{T}^{(\gamma)}\_{\mathfrak{a}}(t) := \mathcal{D}^{\gamma}\_{\star} \mathbf{T}\_{\mathfrak{a}}(t) = [0 \quad \mathcal{D}^{\gamma}\_{\star} \, t^{a} \quad \dots \quad \mathcal{D}^{\gamma}\_{\star} \, t^{aN}].$$

To write the fractional derivative <sup>D</sup>*γ*‹ involved in (1) in the matrix form, the collocation points (15) will be inserted into (18) to have

$$\mathcal{D}^{\gamma}\_{\star} \mathbf{x}\_{N,a}(t\_i) = \mathbf{T}^{(\gamma)}\_a(t\_i) \, \mathbf{D}^t \, \mathbf{A}, \quad i = 0, 1, \ldots, N\_{\prime}$$

which can be expressed equivalently as

$$\mathbf{X}^{(\gamma)} = \mathbf{T}^{(\gamma)} \mathbf{D}^t \mathbf{A},\tag{19}$$

.

where

$$\mathbf{X}^{(\gamma)} = \begin{bmatrix} \mathcal{D}\_{\star}^{\gamma} \mathbf{x}\_{N,a}(t\_0) \\ \mathcal{D}\_{\star}^{\gamma} \mathbf{x}\_{N,a}(t\_1) \\ \vdots \\ \mathcal{D}\_{\star}^{\gamma} \mathbf{x}\_{N,a}(t\_N) \end{bmatrix}, \quad \mathbf{T}^{(\gamma)} = \begin{bmatrix} \mathbf{T}\_{a}^{(\gamma)}(t\_0) \\ \mathbf{T}\_{a}^{(\gamma)}(t\_1) \\ \vdots \\ \mathbf{T}\_{a}^{(\gamma)}(t\_N) \end{bmatrix}.$$

Similarly, the fractional derivative operators D*β<sup>j</sup>* ‹ *x*p*t*q in (1) for *j* " 1, ... , can be approximated as

$$\mathbf{X}^{(\beta\_j)} = \mathbf{T}^{(\beta\_j)} \mathbf{D}^t \mathbf{A},\tag{20}$$

where *X*p*βj*q and *T*p*βj*q are obtained as in (20) by replacing *γ* with *βj*.

By inserting the collocation points into (1), we have the system

$$\mathcal{D}^{\gamma}\_{\star} \mathbf{x}(t\_{i}) = \mathcal{F} \left( t\_{i}, \mathbf{x}(t\_{i}), \mathcal{D}^{\beta\_{1}}\_{\star} \mathbf{x}(t\_{i}), \dots, \mathcal{D}^{\beta\_{\ell}}\_{\star} \mathbf{x}(t\_{i}) \right), \quad i = 0, 1, \dots, N. \tag{21}$$

Considering these equations in a matrix form and substituting the relations (17), (19), and (20) into the resulting system, a fundamental matrix equation is obtained to be solved. Let us assume that the function *F* in (21) is the linear form

$$F = \sum\_{k=1}^{\ell} c\_k(t) \mathcal{D}\_\star^{\mathcal{S}\_k} \mathbf{x}(t) + c\_0(t) \mathbf{x}(t) + h(t),$$

where *ck*p*t*q for *k* " 1, ... , and *<sup>c</sup>*0p*t*q, *h*p*t*q are given functions. In this case, the equations in (21) can be rewritten in the matrix representation as

$$\mathbf{X}^{(\gamma)} = \sum\_{k=1}^{\ell} \mathbf{C}\_k \mathbf{X}^{(\beta\_k)} + \mathbf{C}\_0 \mathbf{X} + \mathbf{H},\tag{22}$$

where the coefficient matrices *Ck*, *k* " 0, 1, ... with size p*N* ` 1qˆp*N* ` 1q and the vector *H* of size p*N* ` 1q ˆ 1 have the forms

$$\mathbf{C}\_{k} = \begin{bmatrix} c\_{k}(t\_{0}) & 0 & \dots & 0 \\ 0 & c\_{k}(t\_{1}) & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & c\_{k}(t\_{N}) \end{bmatrix}, \quad \mathbf{H} = \begin{bmatrix} h(t\_{0}) \\ h(t\_{1}) \\ \vdots \\ h(t\_{N}) \end{bmatrix}'$$

Substituting the relations (17), (19), and (20) into (22), the fundamental matrix equation is obtained

$$\mathbf{W} \mathbf{A} = \mathbf{H}, \quad \text{or} \quad [\mathbf{W}; \mathbf{H}], \tag{23}$$

where

$$\mathbf{W} := \left(\mathbf{T}^{(\gamma)} - \mathbf{C}\_0 \, \mathbf{T} - \mathbf{C}\_1 \, \mathbf{T}^{(\beta\_1)} - \dots - \mathbf{C}\_\ell \, \mathbf{T}^{(\beta\_\ell)}\right) \mathbf{D}^t.$$

Obviously, (23) is a linear matrix equation and *an*, *n* " 0, 1, ... , *N* are the unknowns Bessel coefficients to be determined.

Next aim is to take into account the initial or boundary conditions (2). For the first condition *x*p0q " *x*0, we tend *t* Ñ 0 in (16) to ge<sup>t</sup> the following matrix representation

$$
\hat{\mathbf{X}}\_0 \mathbf{A} = \mathbf{x}\_0, \qquad \hat{\mathbf{X}}\_0 := \mathbf{T}\_a(0) \mathbf{D}^t = \begin{bmatrix} \pounds\_{00} & \pounds\_{01} & \dots & \pounds\_{0N} \end{bmatrix}.
$$

For the remaining initial conditions, one needs to calculate the integer-order derivatives *dk dt<sup>k</sup> <sup>T</sup>α*p*t*q, *k* " 1, 2, ... , *n* ´ 1, which strictly depend on *α* as well as *N*. For example, by choosing *α* " 1{2 and *N* " 7 we ge<sup>t</sup>

$$\mathbf{T}\_{\frac{1}{2}}(t) = \begin{bmatrix} 1 & t^{1/2} & t & t^{3/2} & t^2 & t^{5/2} & t^3 & t^{7/2} \end{bmatrix}.$$

Differentiation twice with respect to *t* reveals that

$$\begin{aligned} \frac{d}{dt}\mathbf{T}\_{\frac{1}{2}}(t) &= \begin{bmatrix} 0 & 0 & 1 & \frac{3}{2}t^{1/2} & 2t & \frac{5}{2}t^{3/2} & 3t^2 & \frac{7}{2}t^{5/2} \end{bmatrix}, \\\ \frac{d^2}{dt^2}\mathbf{T}\_{\frac{1}{2}}(t) &= \begin{bmatrix} 0 & 0 & 0 & 0 & 2 & \frac{15}{4}t^{1/2} & 6t & \frac{35}{4}t^{3/2} \end{bmatrix}. \end{aligned}$$

Now, by differentiating *k* times in (16), and defining

$$\mathbf{T}\_a^{(k)}(t) := \frac{d^k}{dt^k} \mathbf{T}\_a(t),$$

with the limit *t* Ñ 0, we conclude for *k* " 1, 2, . . . , *n* ´ 1 that

$$
\hat{\mathbf{X}}\_k \mathbf{A} = \mathbf{x}\_{k\prime} \qquad \hat{\mathbf{X}}\_k := \mathbf{T}\_a^{(k)}(0) \mathbf{D}^t = [\hat{\mathbf{x}}\_{k0} \quad \hat{\mathbf{x}}\_{k1} \quad \dots \quad \hat{\mathbf{x}}\_{kN}].
$$

Similarly, for the end conditions *x*p*k*qp*L*q " *xLk*, *k* " 0, ... , *n* ´ 1, the following matrix expressions are obtained

$$
\hat{\mathbf{X}}\_{Lk}\mathbf{A} = \mathbf{x}\_{Tk\prime} \qquad \hat{\mathbf{X}}\_{Lk} := \mathbf{T}\_{\mathfrak{a}}^{(k)}(L)\mathbf{D}^{t} = \begin{bmatrix} \pounds\_{L0} & \pounds\_{L1} & \dots & \pounds\_{LN} \end{bmatrix}.
$$

Now, we replace the first *n* rows of the augmented matrix r*W*; *H*s in (23) by the row matrices r*X* p *k*; *xk*s or r*X* p *Lk*; *xLk*s, *k* " 0, 1, . . . , *n* ´ 1 to ge<sup>t</sup> the (nonlinear) algebraic system of equations

$$
\hat{\mathbf{W}}\mathbf{A} = \hat{\mathbf{H}}, \quad \text{or} \quad [\hat{\mathbf{W}}; \hat{\mathbf{H}}].
$$

Thus, the unknown Bessel coefficients in (16) will be known through solving this (nonlinear) system. This can be obtained by using the Newton's iterative algorithm.

**Remark 1.** *In numerical applications below, we frequently encounter the nonlinear terms like x<sup>s</sup>*p*t*q *for s* " 2, 3 ...*. To approximate the nonlinear term x*<sup>2</sup>p*t*q *in terms of <sup>x</sup>*2*N*,*<sup>α</sup>*p*t*q*, the collocation points* (15) *will be substituted into <sup>x</sup>*2*N*,*<sup>α</sup>*p*t*q*. It can be easily seen that in the matrix representation we have*

$$\mathbf{X}^{2} = \begin{bmatrix} \mathbf{x}\_{N,\mathfrak{a}}^{2}(t\_{0}) \\ \mathbf{x}\_{N,\mathfrak{a}}^{2}(t\_{1}) \\ \vdots \\ \mathbf{x}\_{N,\mathfrak{a}}^{2}(t\_{N}) \end{bmatrix} = \begin{bmatrix} \mathbf{x}\_{N,\mathfrak{a}}(t\_{0}) & 0 & \dots & 0 \\ 0 & \mathbf{x}\_{N,\mathfrak{a}}(t\_{1}) & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \mathbf{x}\_{N,\mathfrak{a}}(t\_{N}) \end{bmatrix} \begin{bmatrix} \mathbf{x}\_{N,\mathfrak{a}}(t\_{0}) \\ \mathbf{x}\_{N,\mathfrak{a}}(t\_{1}) \\ \vdots \\ \mathbf{x}\_{N,\mathfrak{a}}(t\_{N}) \end{bmatrix} = \widehat{\mathbf{X}} \mathbf{X}.$$

*Using* (16)*, we further express the matrix X as a product of three block diagonal matrices as follows* p

$$
\hat{\mathbf{X}} = \hat{T}\hat{D}\hat{\mathbf{A}}\_{\prime}
$$

*where*

$$
\hat{\mathbf{T}} = \begin{bmatrix} \mathbf{T}\_a(t\_0) & 0 & \dots & 0 \\ 0 & \mathbf{T}\_a(t\_1) & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \mathbf{T}\_a(t\_N) \end{bmatrix}, \hat{\mathbf{D}} = \begin{bmatrix} \mathbf{D}^t & 0 & \dots & 0 \\ 0 & \mathbf{D}^t & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \mathbf{D}^t \end{bmatrix}, \hat{\mathbf{A}} = \begin{bmatrix} \mathbf{A} & 0 & \dots & 0 \\ 0 & \mathbf{A} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \mathbf{A} \end{bmatrix}.
$$

*Analogously, the higher-order nonlinear terms can be treated recursively Xs* " p*X* p q*<sup>s</sup>*´<sup>1</sup> *X, s* " 3, 4, . . .*.*

#### *4.1. Error Estimation*

In general, the exact solution of most MOFDEs cannot be explicitly obtained. Thus, we need some measurements to test the accuracy of the proposed scheme. Since the truncated Bessel series (11) as an approximate solution is satisfied in (1), our expectation is that the residual error function denoted by <sup>R</sup>*N*,*<sup>α</sup>*p*t*q becomes approximately small. Here, <sup>R</sup>*N*,*<sup>α</sup>*p*t*q : r0, *L*s Ñ R obtained by inserting the computed approximated solution *xN*,*<sup>α</sup>*p*t*q into the differential equation (1). More precisely, for testing accuracy of some numerical models we calculate

$$\mathcal{R}\mathcal{R}\_{\mathbf{N},\mathbf{a}}(t) = \left| \mathcal{D}\_{\star}^{\gamma} \mathbf{x}\_{\mathbf{N},\mathbf{a}}(t) - \mathbb{F} \left( t, \mathbf{x}\_{\mathbf{N},\mathbf{a}}(t), \mathcal{D}\_{\star}^{\beta\_1} \mathbf{x}\_{\mathbf{N},\mathbf{a}}(t), \dots, \mathcal{D}\_{\star}^{\beta\_\ell} \mathbf{x}\_{\mathbf{N},\mathbf{a}}(t) \right) \right| \simeq 0, \quad t \in [0, L]. \tag{24}$$

It should be noticed that the fractional derivatives of order *γ*, *βj*, *j* " 1, ... , of the approximate solution *xN*,*<sup>α</sup>*p*t*q in (24) are calculated by using the properties (5) and (6). Obviously, the residual function is vanished at the collocation points (15), so our expectation is that <sup>R</sup>*N*,*<sup>α</sup>*p*t*q Ñ 0 as *N* tends to infinity. This implies that the smallness of the residual error function shows the closeness of the approximate solution to the true exact solution.

#### **5. Illustrative Test Problems**

Now, we show the benefits of the presented Bessel-collocation scheme by simulating some case examples including various linear and nonlinear initial and boundary value problems. The numerical models and calculations are verified through a comparison with existing computational schemes and experimental measurements. Our computations were carried out using MATLAB software version R2017a.

**Problem 1.** *In the first problem, we consider the following inhomogeneous Bagley–Torvik equation modelling the motion of an immersed plate in a Newtonian fluid [5–7]*

$$\mathbf{x}^{(2)}(t) + \mathcal{D}\_{\bullet}^{\frac{1}{2}}\mathbf{x}(t) + \mathbf{x}(t) = t + 1,$$

*with the initial conditions x*p0q " 1 *and x*p1qp0q " 1*. The exact solution of this problem is x*p*t*q " *t* ` 1*.*

By employing *N* " 2 and *L* " 1, we are looking for an approximate solution in the form *xN*,*<sup>α</sup>*p*t*q " ř<sup>2</sup>*n*"<sup>0</sup> *an*B*<sup>α</sup>n*p*t*q. To this end, we calculate the unknown coefficients *a*0, *a*1, and *a*2. For this example we set *α* " 1 and the collocation points *t*0 " 0, *t*1 " 12 , *t*3 " 0 are used. Using *γ* " 2 and *β*1 " 32, the corresponding matrices and vectors in the fundamental matrix Equation (23) become

$$\mathbf{T}^{(\frac{3}{2})} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 1358/851 \\ 0 & 0 & 167/74 \end{bmatrix}, \quad \mathbf{T}^{(2)} = \begin{bmatrix} 0 & 0 & 2 \\ 0 & 0 & 2 \\ 0 & 0 & 2 \end{bmatrix}, \quad \mathbf{T} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1/2 & 1/4 \\ 1 & 1 & 1 \end{bmatrix}.$$

$$\mathbf{D} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 3 & 3 \end{bmatrix}, \quad \mathbf{H} = \begin{bmatrix} 1 \\ 3/2 \\ 2 \end{bmatrix}, \quad \begin{bmatrix} \hat{\mathbf{W}}; \hat{\mathbf{H}} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \;; & 1 \\ 1 & 3/2 & 1881/134 & \vdots & 3/2 \\ 0 & 1 & 3 & \ddots & 1 \end{bmatrix}.$$

By solving the linear system *W A* p *A* " *H* p , the coefficients matrix is found as

$$\mathbf{A} = \begin{bmatrix} 0 & 1 & 0 \end{bmatrix}^t.$$

Afterwards, by inserting the obtained coefficients into *<sup>x</sup>*2,1p*t*q we ge<sup>t</sup> the approximate solution

$$\mathbf{x}\_{2,1}(t) = \begin{bmatrix} 1 & 1+t & 3t^2+3t+1 \end{bmatrix} \mathbf{A} = 1+t\_{\prime\prime}$$

which is the desired exact solution.

**Problem 2.** *In the next example, the following nonlinear initial-value problem will be considered [5,7]*

$$\mathbf{x}^{(3)}(t) + \mathcal{D}\_{\star}^{\frac{5}{2}}\mathbf{x}(t) + \mathbf{x}^{2}(t) = t^{4}.$$

*The initial conditions are x*p0q " 0*, x*p1qp0q " 0*, and x*p2qp0q " 2*. It can be easily checked that the exact true solution is x*p*t*q " *t*2*.*

For this example, we take *N* " 3, *α* " 1, and the collocation points are t0, 13 , 23 , <sup>1</sup>u. To obtain the unknown coefficients *a*0, *a*1, *a*2, *a*3 in *<sup>x</sup>*3,1p*t*q, the following nonlinear algebraic system of equations to be solved

$$\begin{cases} a\_0 + a\_1 + a\_2 + a\_3 = 0, \\ a\_1 + 3a\_2 + 6a\_3 = 0, \\ 6a\_2 + 30a\_3 = 2, \\ 90a\_3 + \frac{180}{\sqrt{\pi}}a\_3 + (a\_0 + 2a\_1 + 7a\_2 + 37a\_3)^2 = 1. \end{cases}$$

By solving the above system, we ge<sup>t</sup>

$$a\_0 = \frac{2}{3'} \quad a\_1 = -1, \quad a\_2 = \frac{1}{3'} \quad a\_3 = 0.$$

Therefore, we ge<sup>t</sup>

$$\mathbf{x}\_{3,1}(t) = \begin{bmatrix} 1 & 1+t & 1+3t+3t^2 & 1+6t+15t^2+15t^3 \end{bmatrix} \mathbf{A} = t^2, \mathbf{A}$$

which is obviously the exact solution.

In the next example, we show the advantage of using the fractional-order Bessel functions in the computations.

**Problem 3.** *In this test example, we solve the initial-value problem [7]*

$$\mathbf{x}^{(1)}(t) + \mathcal{D}\_{\star}^{\frac{1}{2}}\mathbf{x}(t) + \mathbf{x}(t) = t^{\frac{5}{2}} + \frac{5}{2}t^{\frac{3}{2}} + \frac{15\sqrt{\pi}}{16}t^2,$$

*with initial condition x*p0q " 0*. The exact solution is x*p*t*q " *t*<sup>2</sup>?*t.* *Symmetry* **2020**, *12*, 1260

We first consider *N* " 5 and *α* " 1{2. The approximated solution *<sup>x</sup>*5, 1 2 p*t*q for *t* P r0, 1s takes the form

$$\begin{split} \mathbf{x}\_{5, \frac{1}{2}}(t) &= 3.55309 \times 10^{-14} \, t^2 - 2.20767 \times 10^{-13} \, t + 2.38332 \times 10^{-13} \, t^{1/2} \\ &+ 4.17169 \times 10^{-14} \, t^{3/2} + 1.0 \, t^{5/2} + 9.97959 \times 10^{-111} . \end{split}$$

However, with a lower number of basis functions one can also obtain an accurate result. Using *N* " 2, *α* " 5{2 and *N* " 3, *α* " 5{6, the following approximations are obtained

$$\begin{aligned} \text{tr}\_{2,\frac{5}{7}}(t) &= 4.88118 \times 10^{-17} \, t^5 + 1.0 \, t^{5/2}, \\ \text{tr}\_{3,\frac{5}{7}}(t) &= 1.0 \, t^{5/2} - 4.899133356 \times 10^{-17} \, t^{5/3} + 5.889017161 \times 10^{-17} \, t^{5/6}. \end{aligned}$$

Moreover, to show the advantage of the presented approach and to validate our obtained approximated solutions, we make a comparison in terms of errors in the *L*8 and *L*2 norms in Table 1. We compare the Bessel-collocation approach and the Chelyshkov collocation spectral method [7]. In this comparison, we use different *N* " 1, 2, 3 and *α* " 1{2, 5{2, 5{6.

**Table 1.** Comparison of *L*<sup>8</sup>, *L*2 error norms for test Problem 3.


**Problem 4.** *Consider the boundary value problem [22,23]*

$$
\mathcal{D}^{\gamma}\_{\bullet} \mathbf{x}(t) - \mathcal{D}^{\beta}\_{\bullet} \mathbf{x}(t) = -(1 + \exp \left(t - 1\right)), \quad 1 < \gamma \ll 2, \quad 0 < \beta \ll 1,
$$

*with initial conditions x*p0q " 0 *and x*p1q " 0*. The exact solution corresponds to γ* " 2 *and β* " 1 *is given as x*p*t*q " *t* ´ *t* exp p*t* ´ 1q*.*

Let *N* " 8 and set *α* " 1. For *γ* " 2, *β* " 1, the approximate solution *<sup>x</sup>*8,1p*t*q of the model Problem 4 using Bessel functions in the interval 0 ď *t* ď 1 is

$$\begin{split} \mathbf{x}\_{8,1}(t) &= -0.0001286702494 \, ^\circ \mathbf{f}^{8} - 0.0003938666636 \, ^\circ \mathbf{f}^{\prime} - 0.003196278513 \, ^\circ \mathbf{f}^{6} \\ &- 0.01524130813 \, ^\circ \mathbf{f}^{5} - 0.06134909192 \, \mathbf{f}^{4} - 0.183930672 \, \mathbf{f}^{3} \\ &- 0.3678807668 \, ^\circ \mathbf{f}^{2} + 0.6321206543 \, \mathbf{f} - 2.12897992 \times 10^{-109} \, \text{J} \end{split}$$

In Table 2, we report the numerical results corresponding to these values of *γ*, *β* using different *N* " 8, 16 evaluated at some points *t* P r0, 1s. The corresponding absolute errors E *<sup>N</sup>*,*<sup>α</sup>*p*t*q :" |*x*p*t*q ´ *xN*,*<sup>α</sup>*p*t*q| are also reported in this table. Moreover, the numerical results based on Haar wavelet operational matrices [22] are given in the last column of Table 2. As can see from Table 2, our approximate solutions agree with the results obtained in [22]. The next observation is that more accurate solutions are obtained if one increases the number of Bessel functions *N*.


**Table 2.** Comparison of approximate solutions and the corresponding absolute errors for test Problem 4 using *N* " 8, 16, and *γ* " 2, *β* " 1.

In Figure 1, *<sup>x</sup>*10,1p*t*q is plotted when *γ* " 2 (*β* " 1) is fixed and different values of *β* " 0.25, 0.5, 0.75, 1 (*γ* " 1.25, 1.5, 1.75, 2) are examined. It is observed that as *γ* and *β* approached to 1 and 2 respectively, numerical solutions tend to the exact solutions.

**Figure 1.** Numerical approximations for fixed *β* " 1 and *γ* " 1.25, 1.5, 1.75, 2 (**left**) and fixed *γ* " 2 and *β* " 0.25, 0.5, 0.75, 1 (**right**) in test Problem 4 with *N* " 10.

**Problem 5.** *Let us consider the initial-value problem of Bagley–Torvik equation of fractional order with variable coefficients [24,25]*

$$\mathbf{x}^{(2)}(t) + \frac{1}{2}\sqrt{\pi}t^2\mathcal{D}\_\star^{\frac{3}{2}}\mathbf{x}(t) - 4\sqrt{t}\,\mathbf{x}(t) = 6t\_\star$$

*with initial conditions x*p0q " 0*, x*p1qp0q " 0*. The exact solution is x*p*t*q " *t*3*.*

Clearly, the exact solution is a third-degree polynomial. Therefore, we take *N* " 3 and *α* " 1, which are sufficient to ge<sup>t</sup> the desired approximations. Using the usual collocation points as in Problem 2 and similar to Problem 1, we ge<sup>t</sup> the final augmented matrix

$$
\begin{bmatrix}
\hat{\mathbf{W}}\mathbf{:}\hat{\mathbf{H}}
\end{bmatrix} = \begin{bmatrix}
1 & 1 & 1 & 1 & \vdots & 0 \\
0 & 1 & 3 & 6 & \vdots & 0
\end{bmatrix}.
$$

.

Solving the resulting linear system, we find

$$\mathbf{A} = \begin{bmatrix} -\frac{1}{3} & \frac{3}{5} & -\frac{1}{3} & \frac{1}{15} \end{bmatrix}^t.$$

Therefore, the approximated solution *<sup>x</sup>*3,1p*t*q is obtained as

$$\mathbf{x}\_{3,1}(t) = \begin{bmatrix} 1 & 1+t & 3t^2+3t+1 & 1+6t+15t^2+15t^3 \end{bmatrix} \mathbf{A} = t^3, \mathbf{x}\_{3,1}(t)$$

which is the exact solution.

**Problem 6.** *Consider the fractional Riccati equation [23,26]*

$$\mathcal{D}^{\gamma}\_{\star}\mathfrak{x}(t) + \mathfrak{x}(t) - \mathfrak{x}^{2}(t) = 0, \quad 0 \prec \gamma \lessapprox 1,$$

*on a long time interval with L* " 5 *and initial condition x*p0q " <sup>1</sup>{<sup>2</sup>*. When γ* " 1*, the exact solution is x*p*t*q " 1 expp*t*q`<sup>1</sup> *.*

We calculate the approximated solution *xN*,*<sup>α</sup>*p*t*q using *N* " 7 and *γ* equals to *α* " 1{4. Thus, we ge<sup>t</sup>

$$\begin{split} \mathbf{x}\_{7,\frac{1}{4}}(t) &= 0.1617950181136742 \, t^{\frac{3}{4}} - 0.03445544072182753 \, t^{\frac{1}{2}} - 0.2700823207417999 \, t^{\frac{1}{4}} \\ &- 0.009800299427063008 \, t^{\frac{3}{2}} - 0.1204495580962043 \, t + 0.04675854257899483 \, t^{\frac{5}{4}} \\ &+ 0.0008798927221707359 \, t^{\frac{7}{4}} + 0.49999999999998357401. \end{split}$$

To validate this solution, we also employ the old fractional-order Bessel polynomials as well as Chelyshkov and Legendre functions from the previous works [26,27] with the same parameters as above. The corresponding solutions take the forms respectively

$$\begin{aligned} \mathbf{x}\_{7,\frac{1}{4}}^{B}(t) &= 0.1617932518503192 \, t^{\frac{2}{4}} - 0.03445464899775196 \, t^{\frac{1}{4}} - 0.27008246876491873 \, t^{\frac{5}{4}} \\ &- 0.0097998224671427077 \, t^{\frac{3}{4}} - 0.1204474711992175 \, t + 0.04675716925988139 \, t^{\frac{5}{4}} \\ &+ 0.00087982440038136711651 \, t^{\frac{7}{4}} + 0.5, \end{aligned}$$

$$\begin{aligned} \mathbf{x}\_{7,\frac{1}{4}}^{\mathcal{C}}(t) &= 0.16176395176134591 \, t^{\frac{2}{4}} - 0.034436828673633131 \, t^{\frac{1}{2}} - 0.27008706550355819 \, t^{\frac{1}{4}} \\ &- 0.0097962860256647 \, t^{\frac{2}{4}} - 0.12042146216589336 \, t + 0.046744073064058187 \, t^{\frac{5}{4}} \\ &+ 0.00087942518183550405624 \, t^{\frac{7}{4}} + 0.5, \end{aligned}$$

$$\mathbf{x}\_{7,\frac{1}{4}}^L(t) = 0.16179490530760574 \, t^{\frac{3}{4}} - 0.034455365809483707 \, t^{\frac{1}{2}} - 0.2700823428770175 \, t^{\frac{5}{4}}$$

$$-0.009800287282136090 \, t^{\frac{3}{2}} - 0.12044946374788617 \, t + 0.04675849679405389 \, t^{\frac{5}{4}}$$

$$+ 0.00087989135193320893222 \, t^{\frac{7}{4}} + 0.49999999947316410565$$

To further compare these collocation schemes based on various polynomials, we calculate the estimated residual errors obtained by the relation (24). The graphs of <sup>R</sup>*N*,*<sup>α</sup>*p*t*q on the interval r0, 5s correspond to *γ*, *α* " 1{4 and for *N* " 7 are shown in Figure 2. With respect to Figure 2, it is obviously seen that the residual error functions obtained by the presented Bessel-collocation method are smaller compared to the errors of other polynomial-based numerical collocation schemes.

**Figure 2.** Comparing the error functions in test Problem 6 using old and new Bessel, Chelyshkov, and Legendre functions with *γ*, *α* " 1{4, and *N* " 7.

**Problem 7.** *Consider the following nonlinear boundary value problem with variable coefficients [6]*

$$\ln \mathbf{x}^{(2)}(t) + \Gamma(\frac{4}{5}) \sqrt[5]{t} \mathcal{D}\_{\star}^{\frac{6}{5}} \mathbf{x}(t) + \frac{11}{9} \Gamma(\frac{5}{6}) \sqrt[6]{t} \mathcal{D}\_{\star}^{\frac{1}{6}} \mathbf{x}(t) - \left[ \mathbf{x}^{(1)}(t) \right]^2 = 2 + \frac{1}{10} t^2, \quad 0 < t < 1/2$$

*with boundary conditions x*p0q " 1 *and x*p1q " 2*. The exact solution of this example is x*p*t*q " 1 ` *t*2*.*

In this example, we have *γ* " 2, *β*1 " 6{5, and *β*2 " 1{6. First, we set *α* " 1. The approximate solutions *xN*,*<sup>α</sup>*p*t*q of Problem 7 for *N* " 2, 3 on 0 ď *t* ď 1 are obtained as follows, respectively

$$\begin{aligned} \mathbf{x}\_{2,1}(t) &= 1.0000000000119503851 \, t^2 - 2.98118 \times 10^{-11} \, t + 1.0000000000072342687, \\ \mathbf{x}\_{3,1}(t) &= 1.12820 \times 10^{-9} \, t^3 + 1.0000000012230527969 \, t^2 - 9.66432 \times 10^{-11} \, t \\ &+ 1.0000000000012118324. \end{aligned}$$

To show the gain of the proposed scheme, we compare our results with the collocation method based on Bernstein operational matrix (BOM) of fractional derivative from [6]. Table 3 reports the errors in *L*8 and *L*2 norms of the new Bessel-collocation procedure and the errors of the BOM algorithm. This comparison shows the thoroughness of the proposed method.

**Table 3.** Comparison of *L*<sup>8</sup>, *L*2 error norms for test Problem 7.


**Problem 8.** *We consider the following initial-value problem of multi-term nonlinear fractional differential equation [6]*

$$\mathcal{D}\_{\bullet}^{\gamma} \mathbf{x}(t) + \mathcal{D}\_{\bullet}^{\beta\_1} \mathbf{x}(t) \cdot \mathcal{D}\_{\bullet}^{\beta\_2} \mathbf{x}(t) + [\mathbf{x}(t)]^2 = t^6 + \frac{6t^{3-\gamma}}{\Gamma(4-\gamma)} + \frac{36t^{b-\beta\_1-\beta\_2}}{\Gamma(4-\beta\_1)\Gamma(4-\beta\_2)}, \quad 0 < t < 1, 2$$

*where* 2 ă *γ* ă 3*,* 0 ă *β*1 ă 1*, and* 1 ă *β*2 ă 2 *and the initial conditions are x*p0q " 0*, x*p1qp0q " 0*, and x*p2qp0q " 0*. An easy calculation shows that x*p*t*q " *t*3 *is the exact solution.*

For this example, we set *α* " 1. By applying the collocation technique based upon new Bessel functions at C1: p*γ*, *β*1, *β*2q"p5{2, 9{10, 3{2q and for *N* " 3, 4, the following approximative solutions on 0 ď *t* ď 1 are obtained

$$\begin{split} \mathbf{x}\_{3,1}(t) &= 1.0000000000004519163 \, t^3 - 5.55112 \times 10^{-17} \, t - 5.55112 \times 10^{-17} \, \\ \mathbf{x}\_{4,1}(t) &= 8.60154 \times 10^{-13} \, t^4 + 0.9999999999992085443 \, t^3 - 1.09400 \times 10^{-16} \, t^2 \\ &- 9.83260 \times 10^{-17} \, t + 1.79230 \times 10^{-17} \,. \end{split}$$

A comparison between our collocation scheme at C1 and the method of shifted Jacobi operational matrix (SJOM) [6] with *N* " 24 is made in Table 4. Besides the cases C1 and C2: p*γ*, *β*1, *β*2q " p2.000001, 0.000009, 1.000001q, the following values of p*γ*, *β*1, *β*2q are used in Table 5 for comparison purposes

> C3: p2.99, 0.99, 1.99q, C4: p2.75, 0.75, 1.75q, C5: p2.9999, 0.9999, 1.9999q.

**Table 4.** Comparison of *L*8 error norms for *γ* " 5{2, *β*1 " 9{10, *β*2 " 3{2 in test Problem 8.


**Table 5.** Comparison of *L*8 error norms for various p*γ*, *β*1, *β*2q in test Problem 8.


Looking at Tables 4 and 5 reveals that our numerical solutions obtained via novel Bessel-collocation method are in excellent agreemen<sup>t</sup> with the corresponding exact solutions. Moreover, our proposed scheme is superior compared to the SJOM.

**Problem 9.** *We consider the fractional relaxation-oscillation equation [5,6]*

$$\mathcal{D}\_{\bullet}^{\gamma}x(t) + x(t) = 0, \quad 0 < \gamma < 2,$$

*with the initial condition x*p0q " 1*. If γ* ą 1 *we also have x*p1qp0q " 0*. The exact solution in terms of Mittag–Leffler function is given by x*p*t*q " *<sup>E</sup>γ*p´*t<sup>γ</sup>*q*. Here, <sup>E</sup>γ*p*z*q " ř<sup>8</sup>*k*"<sup>0</sup> *zk* <sup>Γ</sup>p<sup>1</sup>`*γk*q *.*

First, we consider *γ* " 85{100 and set *α* equals to *γ*. We ge<sup>t</sup> the approximated solution *xN*,*<sup>α</sup>*p*t*q using *N* " 8 terms on r0, 1s as follows

$$\begin{cases} \mathbf{x}\_{8,\frac{65}{100}}(t) = 0.0972690897737097 \, t^{\frac{17}{77}} - 0.0264284381134049 \, t^{\frac{17}{4}} + 0.647219778384659 \, t^{\frac{17}{100}} \\ \mathbf{x}\_{1} - 1.05749619232596 \, t^{\frac{17}{100}} + 0.00525953072336809 \, t^{\frac{51}{100}} - 0.284023703239741 \, t^{\frac{51}{100}} \\ \mathbf{x} - 0.000568791172776014 \, t^{\frac{119}{200}} + 1.0 . \end{cases}$$

In Table 6, we calculate the maximum absolute errors using *N* " 8 and *N* " 10. In addition, a comparison is done in this table with the results obtained via SJOM [6]. Looking at Table 6 one can find that the achievement of good approximations to the exact solution is possible using only a few terms of fractional Bessel polynomials.


**Table 6.** Comparison of *L*8 error norms for *γ*, *α* " 85{100 in test Problem 9.

In the next experiments, we investigate the impact of varying *γ* on the maximum absolute errors while *N* " 10 is fixed. Table 7 presents the *L*8 errors for *γ* " 0.2, 0.4, 0.6, 0.8 as well as *γ* " 1.2, 1.4, 1.6, 1.8. In all cases, we exploit *α* " *γ*. Comparisons with existing approximation techniques based on operational matrix of fractional derivatives via B-spline functions [9] and shifted Jacobi functions [6] are also carried out in Table 7.

**Table 7.** Comparison of *L*8 error norms for *N* " 10 and various *γ* in test Problem 9.


**Problem 10.** *In the last case example, let us consider the following singular fractional Lane-Emden type equation [28,29]*

$$\begin{cases} \mathcal{D}\_\star^\gamma \mathbf{x}(t) + \frac{k}{t^{\gamma-\overline{\rho}\_1}} \mathcal{D}\_\star^{\overline{\rho}\_1} \mathbf{x}(t) + \frac{1}{t^{\gamma-2}} \mathbf{x}(t) = \mathbf{g}(t), \quad 0 < t \ll 1, \\\mathbf{x}(0) = 0, \quad \mathbf{x}^{(1)}(0) = 0, \end{cases}$$

*where* 1 ă *γ* ď 2*,* 0 ă *β*1 ď 1*, k* ě 0*, and*

$$\log(t) = t^{2-\gamma} \left( 6t \left( \frac{t^2}{6} + \frac{\Gamma(4-\beta\_1) + k\Gamma(4-\gamma)}{\Gamma(4-\beta\_1)\Gamma(4-\gamma)} \right) - 2 \left( \frac{t^2}{2} + \frac{\Gamma(3-\beta\_1) + k\Gamma(3-\gamma)}{\Gamma(3-\beta\_1)\Gamma(3-\gamma)} \right) \right).$$

*.*

.

*The exact solution is x*p*t*q " *t* 3 ´ *t* 2

To proceed, we take *γ* " 3{2, *β*1 " 1{2, and *k* " 2. Using the collocation points *tj* " 0.001 ` *j*{*N* for *j* " 0, 1, ... , *N* and with *N* " 3, 4, the following approximation solutions are obtained by the Bessel-collocation procedure

$$\begin{aligned} \mathbf{x}\_{3,1}(t) &= 1.0 \, t^3 - 1.0 \, t^2 + 6.23678 \times 10^{-16} \, t + 3.40637 \times 10^{-108} \, , \\ \mathbf{x}\_{4,1}(t) &= -2.93410 \times 10^{-15} \, t^4 + 1.0 \, t^3 - 1.0 \, t^2 + 1.28989 \times 10^{-15} \, t + 3.06104 \times 10^{-108} \, . \end{aligned}$$

Obviously, these approximations are accurate up to machine epsilon. Table 8 reports the comparison of the absolute errors evaluated at some points *t* P r0, 1s obtained by the Bessel-collocation method. For comparison, the results obtained by the collocation method (CM) [29] and the reproducing kernel method (RKM) [28] are also shown in Table 8. The comparisons show that the proposed method is considerably more accurate than other methods.


**Table 8.** Comparison of absolute errors for *γ* " 3{2, *β*1 " 1{2 in test Problem 10.
