Bessel Collocation Method for Solving Fredholm–Volterra Integro-Fractional Differential Equations of Multi-High Order in the Caputo Sense

Abstract

The approximate solutions of Fredholm–Volterra integro-differential equations of multi-fractional order within the Caputo sense (F-VIFDEs) under mixed conditions are presented in this article apply a collocation points technique based completely on Bessel polynomials of the first kind. This new approach depends particularly on transforming the linear equation and conditions into the matrix relations (some time symmetry matrix), which results in resolving a linear algebraic equation with unknown generalized Bessel coefficients. Numerical examples are given to show the technique’s validity and application, and comparisons are made with existing results by applying this process in order to express these solutions, most general programs are written in Python V.3.8.8 (2021).

1. Introduction

Fractional calculus (FC) deals with the differentiation and integration of arbitrary order and it is used in the real world to model and analyze big problems. Fluid flow, electrical networks, fractals theory, control theory, electromagnetic theory, probability, statistics, optics, potential theory, biology, chemistry, diffusion, and viscoelasticity are just a few of the many fields where fractional calculus is used [1,2,3,4].

In recent years, fractional differential equations and integro-fractional differential equations (IFDEs) have captivated the hobby of many researchers in various fields of science and era due to the reality that realistic modeling of a bodily phenomenon with dependencies not only in the immediate time, but also in the past time history can be accomplished effectively using FC. However, in addition to modeling, the solution approaches and their dependability are crucial in detecting key points when a rapid divergence, convergence, or bifurcation begins. As a result, high-precision solutions are always required. Several strategies for solving fractional order differential equations were presented for this purpose (or integro-differential equations), [1,3,4]. The Adomian decomposition method [5], variational iteration method [6], fractional differential transform method [7], fractional difference method [8], and power series method [9] are the most commonly used ideas.

However, from the beginning of 1994, Laguerre, Legendre, Taylor, Fourier, Hermite, and Bessel polynomials have been employed in works [10,11,12,13,14,15] to solve linear differential, integral, and integro-differential difference equations and related systems. In addition, the Bessel polynomial of the first kind method has been used to find approximate solutions of differential, fractional differential equations, integro-differential equations of fractional order, LVIDEs, and LF-VIDEs [16,17,18,19].

The aim of this paper is to expand and apply the first kind of Bessel polynomial in matrix form, as well as the collocation techniques, to evaluate the approximate solution for the multi-high-order linear Fredholm–Volterra integro-fractional differential equations (FVIFDEs) of the general type:

$$\begin{array}{cc}{}_a{}^{c}D_x^{{\sigma }_n}u\left(x\right)\hfill & +{\displaystyle \sum _{l=1}^{n-1}}{\mathit{p}}_l\left(x\right){}_a{}^{c}D_x^{{\sigma }_{n-l}}u\left(x\right)+{\mathit{p}}_n\left(x\right)u\left(x\right) \hfill \\ \hfill & =g\left(x\right)+{\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i{\displaystyle \underset{a}{\overset{b}{\int }}}{F}_i\left(x,t\right){}_a{}^{c}D_t^{{\alpha }_i}u\left(t\right)dt+{\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j{\displaystyle \underset{a}{\overset{x}{\int }}}{V}_j\left(x,t\right){}_a{}^{c}D_t^{{\beta }_j}u\left(t\right)dt, x\in \left[a,b\right]\hfill \end{array}$$

(1)

together with mixed conditions:

$${\displaystyle \sum }_{\ell =0}^{\mu -1}\left\{{𝒽}_{k\ell }{u}^{(\ell )}\left(a\right)+{\overline{𝒽}}_{k\ell }{u}^{(\ell )}\left(b\right)\right\}=C_k, k=0,1,\dots ,\mu -1.$$

(2)

where the fractional orders: ${\sigma }_n>{\sigma }_{n-1}>\cdots >{\sigma }_1>{\sigma }_0=0,{\alpha }_{m_1}>{\alpha }_{m_1-1}>\cdots >{\alpha }_1>{\alpha }_0=0,$ and ${\beta }_{m_2}>{\beta }_{m_2-1}>\cdots >{\beta }_1>{\beta }_0=0,$ and $\mu =\max \left\{\lceil {\sigma }_n\rceil ,\lceil {\alpha }_{m_1}\rceil , \lceil {\beta }_{m_2}\rceil \right\}.$ In addition, $u\left(x\right)$ is an unknown function, the functions ${\mathit{p}}_l\left(x\right), g\left(x\right)\in C\left(\left[a,b\right], ℝ\right), \mathrm{for}\mathrm{all} l=1,2,\dots ,n$, and $F_i\left(x,t\right), V_j\left(x,t\right)\in C\left(S,ℝ\right), \left(\mathrm{with} S=\left\{\left(x,t\right):a\le t\le x\le b\right\}\right)$ are known, with constants ${𝒽}_{k\ell }, {\overline{𝒽}}_{k\ell }, {\lambda }_i,{\overline{\lambda }}_j \mathrm{and} C_k\in ℝ$ for all $k,\ell =0,1,\dots ,\mu -1$, $i=0,1,\dots ,m_1, j=0,1,\dots ,m_2, \left(n, m_1, m_2\in {\mathbb{Z}}^{+}\right)$ are given.

2. Preliminary Considerations

2.1. Basic Definitions and Some Lemmas

Many mathematical definitions of fractional integration and differentiation have come to light in recent years. The most frequently used definitions of fractional calculus involves the Riemann–Liouville fractional derivative and Caputo derivative. In terms of applicability, the Caputo concept is more dependable than the Riemann–Liouville definition. In this section, we are interested some basic definitions and lemmas which are used later on in this paper [1,3,4,20,21].

Definition 1

[22]. A real valued function $u$ defined on closed bounded interval $\left[a,b\right]=I$ be in the space $C_{\gamma }\left(I\right), \gamma \in ℝ,$ if there exist a real number $k>\gamma ,$ such that $u\left(x\right)={\left(x-a\right)}^k{u}_0\left(x\right), where u_0\left(x\right)\in C\left(I\right),$ and it is said to be in the space $C_{\gamma }\left(I\right), \gamma \in ℝ,$ if there exist a real number $k>\gamma ,$ such that $u\left(x\right)={\left(x-a\right)}^k{u}_0\left(x\right), where u_0\left(x\right)\in C\left(I\right),$ and it is said to be in the space $C_{\gamma }^n\left(I\right) iff u^{\left(n\right)}\left(x\right)\in C_{\gamma }\left(I\right), where n\in {\mathbb{Z}}^{+}{{\displaystyle \cup }}^{}\left\{0\right\}.$

Definition 2

[23]. The Riemann–Liouville (R-L) fractional integral operator, ${}_{a}J_x^{\alpha }$, of order $\alpha >0$ of a function $u\in C_{\gamma }\left(I\right),\gamma \ge -1$ is defined as:

$$\begin{array}{cc}{}_{a}J_x^{\alpha }u\left(x\right)\hfill & =\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}\underset{a}{\overset{x}{\int }}{\left(x-t\right)}^{\alpha -1}u\left(t\right)dt, \alpha \in {ℝ}^{+}\hfill \\ {}_{a}J{}_x^{0}u\left(x\right)\hfill & =u\left(x\right).\hfill \end{array}$$

Definition 3

[24]. The Riemann–Liouville (R-L) fractional derivative operator, ${}_a{}^{R}D{}_x^{\alpha }$, of order $\alpha \ge 0$ of a function $u\left(x\right)$ and $u\in C_{-1}^m\left(I\right), m=\alpha$ is normally defined as:

$${}_a{}^{R}D{}_x^{\alpha }u\left(x\right)=D_x^m{}_a{}^{}J{}_x^{m-\alpha }u\left(x\right), m-1 < \alpha \le m, m\in \mathbb{N}.$$

Definition 4

[23]. The Caputo fractional derivative operator, ${}_a{}^{C}D{}_x^{\alpha }$, of a function $u\in C_{-1}^m\left(I\right)$ and $m=\lceil \alpha \rceil$, (ceiling function), is defined as:

$$\begin{array}{cc}{}_a{}^{C}D{}_x^{\alpha }\hfill & ={}_a{}^{}J{}_x^{m-\alpha }\left[D_x^{m}u\right(x\left)\right]\hfill \\ & =\{\begin{array}{cc}\frac{1}{\mathsf{\Gamma }(m-\alpha )}{\int }_a^x{\left(x-t\right)}^{m-\alpha -1}\frac{{\partial }^{m}u\left(t\right)}{\partial t^m}dt,\hfill & m-1<\alpha <m\\ \frac{{\partial }^{m}u\left(x\right)}{\partial x^m},\hfill & \alpha =m,m\in \mathbb{N}\end{array}\hfill \end{array}$$

where the parameter $\alpha$ is the order of the derivative and is allowed to be any positive real number. The operators ${}_{a}J_x^{\alpha }$ and ${}_a{}^{C}D{}_x^{\alpha }$ are linear operators. Furthermore, we have

Lemma 1

[4]. Let $x>a, a\in ℝ$ and for $u\left(x\right)={\left(x-a\right)}^{\beta } for some \beta \ne -m$ is not negative integer, then

$${}_a{}^{ }J{}_x^{\alpha }{\left(x-a\right)}^{\beta }=\frac{\mathsf{\Gamma }\left(\beta +1\right)}{\mathsf{\Gamma }\left(\beta +\alpha +1\right)}{\left(x-a\right)}^{\beta +\alpha }.$$

Lemma 2

[20]. The Caputo derivative of order $\alpha \ge 0$ with $n=\lceil \alpha \rceil$ of the power function $u\left(x\right)={\left(x-a\right)}^{\beta }$ for some $\beta \ge 0$ is formed by:

$${}_a{}^{C}D{}_x^{\alpha }u\left(x\right)=\{\begin{array}{cc}0\hfill & if \beta \in \left\{0,1,2,\cdots ,n-1\right\}\\ \frac{\mathsf{\Gamma }\left(\beta +1\right)}{\mathsf{\Gamma }\left(\beta -\alpha +1\right)}{\left(x-a\right)}^{\beta -\alpha }\hfill & \begin{array}{c}if \beta \in \mathbb{N} and \beta \ge n\\ or \beta \notin \mathbb{N} and \beta >n-1\end{array}\end{array}$$

Lemma 3

[20]. Let $\alpha \ge 0, m=\lceil \alpha \rceil$. Moreover, assume that $u\in C_{-1}^m\left(I\right).$ Then the Caputo fractional derivative ${}_a{}^{C}D{}_x^{\alpha }u\left(x\right)$ is continuous on $I=\left[a,b\right]$ and $\underset{x\to a}{\lim }\left[{}_a{}^{C}D{}_x^{\alpha }u\left(x\right)\right]=0.$

2.2. Bessel Polynomial of the First Kind

The $𝓇$-th degree $N$-truncated Bessel polynomials of the first kind, [25,26], ${\mathrm{J}}_{𝓇}\left(x\right), 𝓇=0,1,\dots ,N$ are defined by

$${\mathrm{J}}_{𝓇}\left(x\right)={\displaystyle \sum }_{\mathrm{𝓀}=0}^{⟦\frac{N-𝓇}{2}⟧}\frac{{\left(-1\right)}^{\mathrm{𝓀}}}{\mathrm{𝓀}!\left(\mathrm{𝓀}+𝓇\right)!}{\left(\frac{x}{2}\right)}^{2\mathrm{𝓀}+𝓇}, 𝓇\in \mathbb{N}, 0\le x<\infty .$$

Here, $N$ is a positive integer that is selected in such a way that $N\ge 𝓇.$ On the other hand, we may express the ${\mathrm{J}}_{𝓇}\left(x\right)$ as follows in the matrix form.

$$J\left(x\right)=\mathit{X}\left(x\right){\mathit{D}}^T or J^T\left(x\right)=\mathit{D}{\mathit{X}}^T\left(x\right).$$

(3)

where $J\left(x\right)=\left[{\mathrm{J}}_0\left(x\right) {\mathrm{J}}_1\left(x\right)\dots {\mathrm{J}}_N\left(x\right)\right]\mathrm{and} \mathit{X}\left(x\right)=\left[1 x x^2 \dots x^N\right]$

If $N$ is odd

$$\mathit{D}={\left[\begin{array}{cccccc}\frac{1}{0!0!2^0}& 0& \frac{-1}{1!1!2^2}& \mathrm{\cdots }& \frac{{(-1)}^{\frac{N-1}{2}}}{\left(\frac{N-1}{2}\right)!\left(\frac{N-1}{2}\right)!2^{N-1}}& 0\\ 0& \frac{1}{0!1!2^1}& 0& \mathrm{\cdots }& 0& \frac{{(-1)}^{\frac{N-1}{2}}}{\left(\frac{N-1}{2}\right)!\left(\frac{N-1}{2}\right)!2^N}\\ 0& 0& \frac{1}{0!2!2^2}& \cdots & \frac{{(-1)}^{\frac{N-3}{2}}}{\left(\frac{N-3}{2}\right)!\left(\frac{N+1}{2}\right)!2^{N-1}}& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & \frac{1}{0!(N-1)!2^{N-1}}& 0\\ 0& 0& 0& \cdots & 0& \frac{1}{0!\mathrm{N}!2^{\mathrm{N}}}\end{array}\right]}_{\left(N+1\right)\times \left(N+1\right)}$$

If $N$ is even

$$\mathit{D}={\left[\begin{array}{cccccc}\frac{1}{0!0!2^0}& 0& \frac{-1}{1!1!2^2}& \mathrm{\cdots }& 0& \frac{{(-1)}^{\frac{N}{2}}}{\left(\frac{N}{2}\right)!\left(\frac{N}{2}\right)!2^N}\\ 0& \frac{1}{0!1!2^1}& 0& \mathrm{\cdots }& \frac{{(-1)}^{\frac{N-2}{2}}}{\left(\frac{N-2}{2}\right)!\left(\frac{N}{2}\right)!2^{N-1}}& 0\\ 0& 0& \frac{1}{0!2!2^2}& \cdots & 0& \frac{{(-1)}^{\frac{N-2}{2}}}{\left(\frac{N+2}{2}\right)!\left(\frac{N+1}{2}\right)!2^N}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & \frac{1}{0!(N-1)!2^{N-1}}& 0\\ 0& 0& 0& \cdots & 0& \frac{1}{0!\mathrm{N}!2^{\mathrm{N}}}\end{array}\right]}_{\left(N+1\right)\times \left(N+1\right)}$$

3. Fundamental Matrix Relations

Recall Equation (1) and rewrite it as follows:

$$D\left({\left\{{\sigma }_l\right\}}_{l=1}^n,x\right)=g\left(x\right)+I_f\left({\left\{{\alpha }_i\right\}}_{i=0}^{m_1},x\right)+I_v\left({\left\{{\beta }_j\right\}}_{j=0}^{m_2},x\right).$$

(4)

where

$$D\left({\left\{{\sigma }_l\right\}}_{l=1}^n,x\right)={}_a{}^{c}D{}_x^{{\sigma }_n}u\left(x\right)+{\displaystyle \sum }_{l=1}^{n-1}{\mathit{p}}_l\left(x\right){}_a{}^{c}D{}_x^{{\sigma }_{n-l}}u\left(x\right)+{\mathit{p}}_n\left(x\right)u\left(x\right)$$

and the integral parts:

$$\begin{array}{cc}{I}_f\left({\left\{{\alpha }_i\right\}}_{i=0}^{m_1},x\right)={\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i\underset{a}{\overset{b}{{\displaystyle \int }}}{F}_i\left(x,t\right){}_a{}^{c}D{}_t^{{\alpha }_i}u\left(t\right)dt,& I_v\left({\left\{{\beta }_j\right\}}_{j=0}^{m_2},x\right)={\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j\underset{a}{\overset{x}{{\displaystyle \int }}}{V}_j\left(x,t\right){}_a{}^{c}D{}_t^{{\beta }_j}u\left(t\right)dt\end{array}.$$

our purpose is to find a close approximation of Equation (1) in the N-truncated Bessel series arrangement

$$u\left(x\right)\cong {\displaystyle \sum }_{𝓇=0}^N{a}_{𝓇}{\mathrm{J}}_{𝓇}\left(x\right).$$

(5)

So that $a_{𝓇}, \mathrm{for}\mathrm{all} 𝓇=0,1,\dots ,N$ are the unknown Bessel coefficients. Before we begin the approximate solution we must convert the solution $u\left(x\right)$ and its ${}_a{}^{c}D{}_x^{{\sigma }_n}u\left(x\right),{}_a{}^{c}D{}_x^{{\sigma }_{n-l}}u\left(x\right), {}_a{}^{c}D{}_x^{{\alpha }_i}u\left(x\right) \mathrm{and} {}_a{}^{c}D{}_x^{{\beta }_j}u\left(x\right), \mathrm{for}\mathrm{all}$ $l=1,2,\dots ,n-1, i=0,1,\dots ,m_1, j=0,1,\dots ,m_2$ in the parts $D\left({\left\{{\sigma }_l\right\}}_{l=1}^n,x\right), I_f\left({\left\{{\alpha }_i\right\}}_{i=0}^{m_1},x\right) \mathrm{and} I_v\left({\left\{{\beta }_j\right\}}_{j=0}^{m_2},x\right),$ to matrix form, within the mixed conditions of Equation (2).

3.1. Matrix Relation for the Fractional Derivative Part $D$

To describe the solution $u\left(x\right)$ of Equation $\left(1\right)$, which is specified by the $N$-truncated Bessel series of Equation $\left(5\right)$. The function defined in relation $\left(5\right)$ in a matrix form

$$\left[u\left(x\right)\right]=J\left(x\right)\mathit{A} ;\mathit{A}={\left[a_0 a_1 a_2\dots a_N\right]}^T$$

(6)

or from Equation $\left(3\right)$

$$\left[u\left(x\right)\right]=\mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}.$$

(7)

The relationship between the matrix $\mathit{X}\left(x\right)$ and its derivative ${\mathit{X}}^{\left(1\right)}\left(x\right)$ is also written as follows:

$${\mathit{X}}^{\left(1\right)}\left(x\right)=\mathit{X}\left(x\right){\mathit{B}}^T.$$

(8)

where

$${\mathit{B}}^T={\left[\begin{array}{cccccc}0& 1& 0& 0& \cdots & 0\\ 0& 0& 2& 0& \cdots & 0\\ 0& 0& 0& 3& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& 0& \cdots & N\\ 0& 0& 0& 0& \cdots & 0\end{array}\right]}_{\left(N+1\right)\times \left(N+1\right)}$$

We will also get the recurrence relations from Equation (8):

$$\begin{gathered} {\mathit{X}}^{\left(0\right)}\left(x\right)=\mathit{X}\left(x\right) \\ {\mathit{X}}^{\left(1\right)}\left(x\right)=\mathit{X}\left(x\right){\mathit{B}}^T \\ {\mathit{X}}^{\left(2\right)}\left(x\right)={\mathit{X}}^{\left(1\right)}\left(x\right){\mathit{B}}^T=\mathit{X}\left(x\right){\left({\mathit{B}}^T\right)}^2 \\ ⋮ \\ {\mathit{X}}^{\left(𝒿\right)}\left(x\right)={\mathit{X}}^{\left(𝒿-1\right)}\left(x\right){\mathit{B}}^T=\mathit{X}\left(x\right){\left({\mathit{B}}^T\right)}^{𝒿}. \end{gathered}$$

(9)

Here, note that ${\left({\mathit{B}}^T\right)}^0={\left[I\right]}_{\left(N+1\right)\times \left(N+1\right)}$ is an identity matrix of dimension $\left(N+1\right)$. Using mathematical induction, we can prove that Equation (9) is correct. By applying the same concept to Equation (7) and using Equation (9), we attain matrix relation

$$\begin{array}{cc}{u}^{\left(𝒿\right)}\left(x\right)\hfill & ={\mathit{X}}^{\left(𝒿\right)}\left(x\right){\mathit{D}}^T\mathit{A}\hfill \\ u^{\left(𝒿\right)}\left(x\right)\hfill & =\mathit{X}\left(x\right){\left({\mathit{B}}^T\right)}^{𝒿}{\mathit{D}}^T\mathit{A},\mathrm{for}\mathrm{each}𝒿=0,1,\dots ,\mu \mathrm{and} \mu =\max \left\{\lceil {\sigma }_n\rceil ,\lceil {\alpha }_{m_1}\rceil , \lceil {\beta }_{m_2}\rceil \right\}\end{array}$$

(10)

By using Equation (7) with (9) and applying the Caputo Definition 4, with Lemma 1 and 2, we can convert the fractional terms ${}_a{}^{c}D{}_x^{{\sigma }_{n-l}}u\left(x\right),\overline{n}\left({\sigma }_{n-l}\right)-1<{\sigma }_{n-l}\le \overline{n}\left({\sigma }_{n-l}\right),$ that is $\overline{n}\left({\sigma }_{n-l}\right)=\lceil {\sigma }_{n-l}\rceil ,$ for all $l=0,1,\dots ,n-1$ to matrix form:

$$\begin{array}{cc}{}_a{}^{c}D{}_x^{{\sigma }_{n-l}}u\left(x\right)\hfill & ={}_a{}^{c}D{}_x^{{\sigma }_{n-l}}\mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={}_a{}^{ }J{}_x^{\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)}{D}^{ \overline{n}\left({\sigma }_{n-l}\right)} \mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={}_a{}^{ }J{}_x^{\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)} \mathit{X}\left(x\right){\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}{\mathit{D}}^T\mathit{A}\hfill \\ \hfill & =x^{\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}} \mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)}{\mathit{D}}^T\mathit{A}.\hfill \end{array}$$

Since

$${}_a{}^{ }J{}_x^{\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)} \mathit{X}\left(x\right)={}_a{}^{ }J{}_x^{\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)}\left[1 x x^2\dots x^N\right]$$

$$\begin{array}{cc}\hfill =[\frac{\mathsf{\Gamma }\left(1\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+1\right)}{x}^{\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+0},& \frac{\mathsf{\Gamma }\left(2\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+2\right)}{x}^{\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+1},\hfill \\ \hfill \dots ,& \frac{\mathsf{\Gamma }\left(\mathrm{N}+1\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+N+1\right)}{x}^{\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+N}]\hfill \end{array}$$

$$=x^{\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}}\left[1x{x}^2\dots x^N\right]\left[\begin{array}{cccc}\frac{\mathsf{\Gamma }\left(1\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+1\right)}& 0& \cdots & 0\\ 0& \frac{\mathsf{\Gamma }\left(2\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+2\right)}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{\mathsf{\Gamma }\left(\mathrm{N}+1\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+N+1\right)}\end{array}\right]$$

Putting

$$\mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)=\left[\begin{array}{cccc}\frac{\mathsf{\Gamma }\left(1\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+1\right)}& 0& \cdots & 0\\ 0& \frac{\mathsf{\Gamma }\left(2\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+2\right)}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{\mathsf{\Gamma }\left(\mathrm{N}+1\right)}{\mathsf{\Gamma }\left( \overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}+N+1\right)}\end{array}\right]$$

(11)

Thus, for all $l=1,2,\dots ,n-1$ in general we obtain

$${}_a{}^{c}D{}_x^{{\sigma }_{n-l}}u\left(x\right)={}_a{}^{c}D{}_x^{{\sigma }_{n-l}}\mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}=x^{\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}} \mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}{\mathit{D}}^T\mathit{A}.$$

(12)

and

$${}_a{}^{c}D{}_x^{{\sigma }_n}u\left(x\right)={}_a{}^{c}D{}_x^{{\sigma }_n}\mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}=x^{\overline{n}\left({\sigma }_n\right)-{\sigma }_n} \mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\sigma }_n\right)-{\sigma }_n\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)}{\mathit{D}}^T\mathit{A}.$$

(13)

Using mathematical induction, we can prove that Equations (12) and (13) are correct. By substituting expressions $\left(7\right), \left(12\right)$ and $\left(13\right)$ into $\left(4\right)$, As well we can make this assumption $\mathit{y}\left(n,x\right)=x^{\overline{n}\left({\sigma }_n\right)-{\sigma }_n}, \mathit{y}\left(n-l,x\right)=x^{\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}},\mathrm{for}\mathrm{all} l=1,2,\dots ,n-1$, we have

$$\begin{array}{cc}D\left({\left\{{\sigma }_l\right\}}_{l=1}^n,x\right)\hfill & =[\mathit{y}\left(n,x\right) \mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\sigma }_n\right)-{\sigma }_n\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)}{\mathit{D}}^T\hfill \\ \hfill & +{\displaystyle \sum _{l=1}^{n-1}}{\mathit{p}}_l\left(x\right)\mathit{y}\left(n-l,x\right) \mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}{\mathit{D}}^T+{\mathit{p}}_n\left(x\right)\mathit{X}\left(x\right){\mathit{D}}^T]\mathit{A}.\hfill \end{array}$$

(14)

3.2. Matrix Relation for the Fredholm Integral Part $I_f$

The $N$-truncated Taylor series around (0,0), [27] and the $N$-truncated Bessel series can be used to approximate the Fredholm kernel functions $F_i\left(x,t\right), i=0,1,\dots ,m_1$, respectively

$$F_i\left(x,t\right)={\displaystyle \sum }_{𝓂=0}^N{\displaystyle \sum }_{𝓃=0}^N{}_t{}^{ }F{}_{𝓂𝓃}^i{x}^{𝓂}{t}^{𝓃} \mathrm{and} F_i\left(x,t\right)={\displaystyle \sum }_{𝓂=0}^N{\displaystyle \sum }_{𝓃=0}^N{}_b{}^{ }F{}_{𝓂𝓃}^i{\mathrm{J}}_{𝓂}\left(x\right){\mathrm{J}}_{𝓃}\left(t\right), i=0,1,\dots ,m_1$$

(15)

where

$$\left[{}_t{}^{ }F{}_{𝓂𝓃}^i\right]=\frac{1}{𝓂!𝓃!}\frac{{\partial }^{𝓂+𝓃}{F}_i\left(0,0\right)}{\partial x^{𝓂}\partial t^{𝓃}}, i=0,1,\dots ,m_1, 𝓂,𝓃=0,1,\dots ,N.$$

In matrix forms, the Equation (15) may be written as Equations (16) and (17), respectively

$$F_i\left(x,t\right)=\mathit{X}\left(x\right){\mathit{F}}_t^i{\mathit{X}}^T\left(t\right), {\mathit{F}}_t^i=\left[{}_t{}^{ }F{}_{𝓂𝓃}^i\right], i=0,1,\dots ,m_1, 𝓂,𝓃=0,1,\dots ,N.$$

(16)

and

$$F_i\left(x,t\right)=J\left(x\right){\mathit{F}}_b^i{J}^T\left(t\right),{\mathit{F}}_b^i=\left[{}_b{}^{ }F{}_{𝓂𝓃}^i\right], i=0,1,\dots ,m_1, 𝓂,𝓃=0,1,\dots ,N.$$

(17)

From Equations (16) and (17), it also comes out according to Equation (3), the following relation

$$\begin{gathered} \mathit{X}\left(x\right){\mathit{F}}_t^i{\mathit{X}}^T\left(t\right)=J\left(x\right){\mathit{F}}_b^i{J}^T\left(t\right), i=0,1,\dots ,m_1 \\ \mathit{X}\left(x\right){\mathit{F}}_t^i{\mathit{X}}^T\left(t\right)=\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{F}}_b^i\mathit{D}{\mathit{X}}^T\left(t\right), i=0,1,\dots ,m_1 \\ {\mathit{F}}_t^i={\mathit{D}}^T{\mathit{F}}_b^i\mathit{D} \mathrm{or} {\mathit{F}}_b^i={\left({\mathit{D}}^T\right)}^{-1}{\mathit{F}}_t^i{\mathit{D}}^{-1}, i=0,1,\dots ,m_1 \end{gathered}$$

(18)

In the same way from Equations (12) and (13), convert ${}_a{}^{c}D{}_x^{{\alpha }_i}u\left(x\right),\mathrm{and} \overline{n}\left({\alpha }_i\right)-1<{\alpha }_i\le \overline{n}\left({\alpha }_i\right), \mathrm{i}.\mathrm{e}., \overline{n}\left({\alpha }_i\right)=\lceil {\alpha }_i\rceil \mathrm{for}\mathrm{all} i=0,1,\dots ,m_1,$ by apply the Caputo Definition 4 with Lemma 1 to the matrix form, we obtain

$$\begin{array}{cc}{}_a{}^{c}D{}_x^{{\alpha }_i}u\left(x\right)\hfill & ={}_a{}^{c}D{}_x^{{\alpha }_i}\mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={}_a{}^{ }J{}_x^{\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right)}{D}^{ \overline{n}\left({\alpha }_i\right)} \mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={}_a{}^{ }J{}_x^{\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right)} \mathit{X}\left(x\right){\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}\hfill \\ \hfill & =x^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i} \mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}, i=0,1,\dots ,m_1\hfill \end{array}$$

(19)

where $\mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right)$ is defined at Equation (11). We obtain the matrix relation (20), put the Equation (3) into (17), and then replace the obtained matrix with matrix (19) in the Fredholm integral part $I_f$ in Equation (4)

$$\begin{array}{cc}{I}_f\left({\left\{{\alpha }_i\right\}}_{i=0}^{m_1},x\right)\hfill & ={\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i{\displaystyle \underset{a}{\overset{b}{\int }}}{F}_i\left(x,t\right){}_a{}^{c}D{}_t^{{\alpha }_i}u\left(t\right)dt\hfill \\ \hfill & ={\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i{\displaystyle \underset{a}{\overset{b}{\int }}}J\left(x\right){\mathit{F}}_b^i{J}^T\left(t\right)t^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i} \mathit{X}\left(t\right) \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}dt\hfill \\ \hfill & ={\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i{\displaystyle \underset{a}{\overset{b}{\int }}}\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{F}}_b^i\mathit{D}{\mathit{X}}^T\left(t\right)t^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i} \mathit{X}\left(t\right) \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}dt\hfill \\ \hfill & ={\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{F}}_b^i\mathit{D}\left({\displaystyle \underset{a}{\overset{b}{\int }}}{X}^T\left(t\right) \mathit{X}\left(t\right)t^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i} dt\right)\mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{F}}_b^i\mathit{D} {\mathit{H}}_{f,i} \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}.\hfill \end{array}$$

(20)

where

$$\begin{gathered} {\mathit{H}}_{f,i}=\underset{a}{\overset{b}{{\displaystyle \int }}}{\mathit{X}}^T\left(t\right) \mathit{X}\left(t\right)t^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i} dt=\left[h_{rs}^{f,i}\right], i=0,1,\dots ,m_1, r,s=0,1,\dots ,N. \\ \left[h_{rs}^{f,i}\right]=\frac{b^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i+r+s+1}-a^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i+r+s+1}}{\overline{n}\left({\alpha }_i\right)-{\alpha }_i+r+s+1}, i=0,1,\dots ,m_1, r,s=0,1,\dots ,N. \end{gathered}$$

We can get the last matrix form (21) by replacing the matrix relation (18) into expression (20).

$$I_f\left({\left\{{\alpha }_i\right\}}_{i=0}^{m_1},x\right)={\displaystyle \sum }_{i=0}^{m_1}{\lambda }_i\mathit{X}\left(x\right){\mathit{F}}_t^i {\mathit{H}}_{f,i} \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}.$$

(21)

3.3. Matrix Relation for the Volterra Integral Part $I_v$

The $N$-truncated Taylor series around (0,0), [27] and the $N$-truncated Bessel series can be used to approximate the Volterra kernel functions $V_j\left(x,t\right), j=0,1,\dots ,m_2$, respectively

$$V_j\left(x,t\right)={\displaystyle \sum }_{𝓂=0}^N{\displaystyle \sum }_{𝓃=0}^N{}_t{}^{ }V{}_{𝓂𝓃}^j{x}^{𝓂}{t}^{𝓃} \mathrm{and} V_j\left(x,t\right)={\displaystyle \sum }_{𝓂=0}^N{\displaystyle \sum }_{𝓃=0}^N{}_b{}^{ }V{}_{𝓂𝓃}^j{\mathrm{J}}_{𝓂}\left(x\right){\mathrm{J}}_{𝓃}\left(t\right), j=0,1,\dots ,m_2$$

(22)

where

$$\left[{}_t{}^{ }V{}_{𝓂𝓃}^j\right]=\frac{1}{𝓂!𝓃!}\frac{{\partial }^{𝓂+𝓃}{V}_j\left(0,0\right)}{\partial x^{𝓂}\partial t^{𝓃}}, j=0,1,\dots ,m_2, 𝓂,𝓃=0,1,\dots ,N.$$

The relations in Equation (22) can be transformed into matrix forms:

$$V_j\left(x,t\right)=\mathit{X}\left(x\right){\mathit{V}}_t^j{\mathit{X}}^T\left(t\right), {\mathit{V}}_t^j=\left[{}_t{}^{ }V{}_{𝓂𝓃}^j\right], j=0,1,\dots ,m_2, 𝓂,𝓃=0,1,\dots ,N$$

(23)

and

$$V_j\left(x,t\right)=J\left(x\right){\mathit{V}}_b^j{J}^T\left(t\right), {\mathit{V}}_b^j=\left[{}_b{}^{ }V{}_{𝓂𝓃}^j\right], j=0,1,\dots ,m_2, 𝓂,𝓃=0,1,\dots ,N$$

(24)

from Equations (23) and (24), it also comes out according to Equation (3) we obtain the following relation:

$$\begin{gathered} \mathit{X}\left(x\right){\mathit{V}}_t^j{\mathit{X}}^T\left(t\right)=J\left(x\right){\mathit{V}}_b^j{J}^T\left(t\right), j=0,1,\dots ,m_2 \\ \mathit{X}\left(x\right){\mathit{V}}_t^j{\mathit{X}}^T\left(t\right)=\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{V}}_b^j\mathit{D}{\mathit{X}}^T\left(t\right), j=0,1,\dots ,m_2 \\ {\mathit{V}}_t^j={\mathit{D}}^T{\mathit{V}}_b^j\mathit{D} \mathrm{or} {\mathit{V}}_b^j={\left({\mathit{D}}^T\right)}^{-1}{\mathit{V}}_t^j{\mathit{D}}^{-1}, j=0,1,\dots ,m_2 \end{gathered}$$

(25)

Finally, in the same way from Equations (12) and (13), convert ${}_a{}^{c}D{}_x^{{\beta }_j}u\left(x\right), \overline{n}\left({\beta }_j\right)-1<{\beta }_j\le \overline{n}\left({\beta }_j\right), i.e., \overline{n}\left({\beta }_j\right)=\lceil {\beta }_j\rceil , \mathrm{for}\mathrm{all} j=0,1,\dots ,m_2$, by applying the Caputo Definition 4 with Lemma 1, 2 and 3 to matrix form, we obtain

$$\begin{array}{cc}{}_a{}^{c}D{}_x^{{\beta }_j}u\left(x\right)\hfill & ={}_a{}^{c}D{}_x^{{\beta }_j}\mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={}_a{}^{ }J{}_x^{\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right)}{D}^{ \overline{n}\left({\beta }_j\right)} \mathit{X}\left(x\right){\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={}_a{}^{ }J{}_x^{\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right)} \mathit{X}\left(x\right){\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}\hfill \\ \hfill & =x^{\overline{n}\left({\beta }_j\right)-{\beta }_j} \mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}, j=0,1,\dots ,m_2\hfill \end{array}$$

(26)

where $\mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right)$ define at Equation (11). We obtain the matrix relation (27), put the Equation (3) into (24) and then replace the obtained matrix with matrix (26) in Fredholm integral part $I_v$ in Equation (4)

$$\begin{array}{cc}{I}_v\left({\left\{{\beta }_j\right\}}_{j=0}^{m_2},x\right)\hfill & ={\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j{\displaystyle \underset{a}{\overset{x}{\int }}}{V}_j\left(x,t\right){}_a{}^{c}D{}_t^{{\beta }_j}u\left(t\right)dt\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j{\displaystyle \underset{a}{\overset{x}{\int }}}J\left(x\right){\mathit{V}}_b^j{J}^T\left(t\right)t^{\overline{n}\left({\beta }_j\right)-{\beta }_j} \mathit{X}\left(t\right) \mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}dt\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j{\displaystyle \underset{a}{\overset{x}{\int }}}\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{V}}_b^j\mathit{D}{\mathit{X}}^T\left(t\right)t^{\overline{n}\left({\beta }_j\right)-{\beta }_j} \mathit{X}\left(t\right) \mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}dt\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{V}}_b^j\mathit{D}\left({\displaystyle \underset{a}{\overset{x}{\int }}}{\mathit{X}}^T\left(t\right) \mathit{X}\left(t\right)t^{\overline{n}\left({\beta }_j\right)-{\beta }_j}dt\right) \mathit{C}(\overline{n}\left({\beta }_j\right)\hfill \\ \hfill & -{\beta }_j) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j\mathit{X}\left(x\right){\mathit{D}}^T{\mathit{V}}_b^j\mathit{D}{\mathit{H}}_{v,j}\left(x\right)\mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}.\hfill \end{array}$$

(27)

where

$$\begin{gathered} {\mathit{H}}_{v,j}\left(x\right)=\underset{a}{\overset{x}{{\displaystyle \mathsf{\int }}}}{\mathit{X}}^T\left(t\right) \mathit{X}\left(t\right)t^{\overline{n}\left({\beta }_j\right)-{\beta }_j}dt=\left[h_{rs}^{v,j}\left(x\right)\right], j=0,1,\dots ,m_2, r,s=0,1,\dots ,N. \\ \left[h_{rs}^{v,j}\left(x\right)\right]=\frac{x^{\overline{n}\left({\beta }_j\right)-{\beta }_j+r+s+1}-a^{\overline{n}\left({\beta }_j\right)-{\beta }_j+r+s+1}}{\overline{n}\left({\beta }_j\right)-{\beta }_j+r+s+1}, j=0,1,\dots ,m_2, r,s=0,1,\dots ,N. \end{gathered}$$

We can get the last matrix form (28) by replacing the matrix relation (25) into expression (27).

$$I_v\left({\left\{{\beta }_j\right\}}_{j=0}^{m_2},x\right)={\displaystyle \sum }_{j=0}^{m_2}{\overline{\lambda }}_j\mathit{X}\left(x\right){\mathit{V}}_t^j{\mathit{H}}_{v,j}\left(x\right)\mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}.$$

(28)

3.4. Matrix Relation for the Conditions

For each $k=0,1,\dots ,\mu -1$ and $\mu =max\left\{\lceil {\sigma }_n\rceil ,\lceil {\alpha }_{m_1}\rceil ,\lceil {\beta }_{m_2}\rceil \right\},$ applying the relation (10) to each mixed condition of Equation (2), we obtain the corresponding condition matrix forms as follows.

$${\displaystyle \sum }_{\ell =0}^{\mu -1}\left\{{𝒽}_{k\ell }\mathit{X}\left(a\right){\left({\mathit{B}}^T\right)}^{\ell }{\mathit{D}}^T\mathit{A}+{\overline{𝒽}}_{k\ell }\mathit{X}\left(b\right){\left({\mathit{B}}^T\right)}^{\ell }{\mathit{D}}^T\mathit{A}\right\}=\left[{\mathit{C}}_k\right], k=0,1,\dots ,\mu -1$$

Thus

$${\displaystyle \sum }_{\ell =0}^{\mu -1}\left[{𝒽}_{k\ell }\mathit{X}\left(a\right)+{\overline{𝒽}}_{k\ell }\mathit{X}\left(b\right)\right]{\left({\mathit{B}}^T\right)}^{\ell }{\mathit{D}}^T\mathit{A}=\left[{\mathit{C}}_k\right], k=0,1,\dots ,\mu -1.$$

(29)

4. Method of Solution

To construct the fundamental matrix equation that corresponds to Equation (1), insert the matrix relations (14), (21), and (28) into Equation (4) to obtain the following matrix equation

$$\begin{gathered} [\mathit{y}\left(n,x\right)\mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\sigma }_n\right)-{\sigma }_n\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)}{\mathit{D}}^T+{\displaystyle \sum _{l=1}^{n-1}}{\mathit{p}}_l\left(x\right)\mathit{y}\left(n-l,x\right)\mathit{X}\left(x\right) \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}{\mathit{D}}^T \\ \begin{array}{cc}\hfill & +{\mathit{p}}_n\left(x\right)\mathit{X}\left(x\right){\mathit{D}}^T]\mathit{A}\hfill \\ \hfill & =g\left(x\right)+{\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i\mathit{X}\left(x\right){\mathit{F}}_t^i {\mathit{H}}_{f,i} \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j\mathit{X}\left(x\right){\mathit{V}}_t^j{\mathit{H}}_{v,j}\left(x\right)\mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}.\hfill \end{array} \end{gathered}$$

(30)

We get the following system of equations by setting the collocation points, [28], described by $x_{\mathrm{𝒾}}=a+\frac{b-a}{N}\mathrm{𝒾}, \mathrm{𝒾}=0,1,\dots ,N$:

$$\begin{gathered} [\mathit{y}\left(n,x_{\mathrm{𝒾}}\right) \mathit{X}\left(x_{\mathrm{𝒾}}\right) \mathit{C}\left(\overline{n}\left({\sigma }_n\right)-{\sigma }_n\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)}{\mathit{D}}^T+{\displaystyle \sum _{l=1}^{n-1}}{\mathit{p}}_l\left(x_{\mathrm{𝒾}}\right)\mathit{y}\left(n-l,x_{\mathrm{𝒾}}\right) \mathit{X}\left(x_{\mathrm{𝒾}}\right) \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}{\mathit{D}}^T \\ \begin{array}{cc}\hfill & +{\mathit{p}}_n\left(x_{\mathrm{𝒾}}\right)\mathit{X}\left(x_{\mathrm{𝒾}}\right){\mathit{D}}^T]\mathit{A}\hfill \\ \hfill & =g\left(x_{\mathrm{𝒾}}\right)+{\displaystyle \sum _{i=0}^{m_1}}{\lambda }_i\mathit{X}\left(x_{\mathrm{𝒾}}\right){\mathit{F}}_t^i {\mathit{H}}_{f,i} \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T\mathit{A}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{m_2}}{\overline{\lambda }}_j\mathit{X}\left(x_{\mathrm{𝒾}}\right){\mathit{V}}_t^j{\mathit{H}}_{v,j}\left(x_{\mathrm{𝒾}}\right)\mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}{\mathit{D}}^T\mathit{A}, \mathrm{𝒾}=0,1,\dots ,N.\hfill \end{array} \end{gathered}$$

or in brief, the most important matrix equation is

$$\begin{gathered} [\mathit{y}\left(n\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_n\right)-{\sigma }_n\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)}{\mathit{D}}^T+{\displaystyle \sum }_{l=1}^{n-1}{\mathit{p}}_l\mathit{y}\left(n-l\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}{\mathit{D}}^T+{\mathit{p}}_n\mathit{X}{\mathit{D}}^T \\ -{\displaystyle \sum }_{i=0}^{m_1}{\lambda }_i\mathit{X}{\mathit{F}}_t^i {\mathit{H}}_{f,i} \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T-{\displaystyle \sum }_{j=0}^{m_2}{\overline{\lambda }}_j\overline{\mathit{X}}{\overline{\mathit{V}}}^j{\overline{\mathit{H}}}_j{\overline{\mathit{C}}}_j{\overline{\mathit{B}}}^j\overline{\mathit{D}}]\mathit{A}=\mathit{G}. \end{gathered}$$

(31)

where

$$\mathit{y}\left(n-l\right)=\left[\begin{array}{cccc}{x}_0^{\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)}& 0& \cdots & 0\\ 0& x_1^{\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & x_N^{\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right)}\end{array}\right], {\mathit{p}}_l=\left[\begin{array}{cccc}{p}_l\left(x_0\right)& 0& \cdots & 0\\ 0& p_l\left(x_1\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & p_l\left(x_N\right)\end{array}\right]$$

for all $l=0,1,\dots ,n,\mathrm{and} \mathit{y}\left(0\right)={\left[I\right]}_{\left(N+1\right)\times \left(N+1\right)}, {\mathit{p}}_0={\left[I\right]}_{\left(N+1\right)\times \left(N+1\right)}$ are the unit matrix,

$$\mathit{X}=\left[\begin{array}{c}\mathit{X}\left(x_0\right)\\ \mathit{X}\left(x_1\right)\\ ⋮\\ \mathit{X}\left(x_N\right)\end{array}\right]=\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^N\\ 1& x_1& x_1^2& \cdots & x_1^N\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ 1& x_N& x_N^2& \cdots & x_N^N\end{array}\right], \overline{\mathit{X}}=\left[\begin{array}{cccc}\mathit{X}\left(x_0\right)& 0& \cdots & 0\\ 0& \mathit{X}\left(x_1\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \mathit{X}\left(x_N\right)\end{array}\right]$$

also, for each fractional order $\gamma =\left\{{\sigma }_n,{\sigma }_{n-l},{\alpha }_i\mathrm{and}{\beta }_j,\mathrm{for}\mathrm{all}l=1,2,\dots ,n,i=0,1,\dots ,m_1,j=0,1,\dots ,m_2\right\}$ we are putting

$$\begin{gathered} \mathit{C}\left(\overline{n}\left(\gamma \right)-\gamma \right)=\left[\begin{array}{cccc}\frac{\mathsf{\Gamma }\left(1\right)}{\mathsf{\Gamma }\left(\overline{n}\left(\gamma \right)-\gamma +1\right)}& 0& \cdots & 0\\ 0& \frac{\mathsf{\Gamma }\left(2\right)}{\mathsf{\Gamma }\left(\overline{n}\left(\gamma \right)-\gamma +2\right)}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{\mathsf{\Gamma }\left(N+1\right)}{\mathsf{\Gamma }\left(\overline{n}\left(\gamma \right)-\gamma +\mathrm{N}+1\right)}\end{array}\right] \\ {\left({\mathit{B}}^T\right)}^{\overline{n}\left(\gamma \right)}={\left[\begin{array}{cc}\begin{array}{ccc}0& 1& 0\\ 0& 0& 2\\ 0& 0& 0\end{array}& \begin{array}{ccc}0& \dots & 0\\ 0& \dots & 0\\ 3& \dots & 0\end{array}\\ \begin{array}{ccc}⋮& ⋮& ⋮\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}⋮& \ddots & ⋮\\ 0& \dots & N\\ 0& \dots & 0\end{array}\end{array}\right]}^{\overline{n}\left(\gamma \right)},{\mathit{F}}_t^i=\left[\begin{array}{cccc}{}_t{}^{ }F{}_{00}^i& {}_t{}^{ }F{}_{01}^i& \cdots & {}_t{}^{ }F{}_{0N}^i\\ {}_t{}^{ }F{}_{10}^i& {}_t{}^{ }F{}_{11}^i& \cdots & {}_t{}^{ }F{}_{1N}^i\\ ⋮& ⋮& \ddots & ⋮\\ {}_t{}^{ }F{}_{N0}^i& {}_t{}^{ }F{}_{N1}^i& \cdots & {}_t{}^{ }F{}_{NN}^i\end{array}\right],{\mathit{V}}_t^j=\left[\begin{array}{cccc}{}_t{}^{ }V{}_{00}^j& {}_t{}^{ }V{}_{01}^j& \cdots & {}_t{}^{ }V{}_{0N}^j\\ {}_t{}^{ }V{}_{10}^j& {}_t{}^{ }V{}_{11}^j& \cdots & {}_t{}^{ }V{}_{1N}^j\\ ⋮& ⋮& \ddots & ⋮\\ {}_t{}^{ }V{}_{N0}^j& {}_t{}^{ }V{}_{N1}^j& \cdots & {}_t{}^{ }V{}_{NN}^j\end{array}\right] \\ i=0,1,\dots ,m_1, j=0,1,\dots ,m_2 \end{gathered}$$

where

$$\begin{gathered} \left[{}_t{}^{ }F{}_{𝓂𝓃}^i\right]=\frac{1}{𝓂!𝓃!}\frac{{\partial }^{𝓂+𝓃}{F}_i\left(0,0\right)}{\partial x^{𝓂}\partial t^{𝓃}}, i=0,1,\dots ,m_1, 𝓂,𝓃=0,1,\dots ,N \\ \left[{}_t{}^{ }V{}_{𝓂𝓃}^j\right]=\frac{1}{𝓂!𝓃!}\frac{{\partial }^{𝓂+𝓃}{V}_j\left(0,0\right)}{\partial x^{𝓂}\partial t^{𝓃}}, j=0,1,\dots ,m_2, 𝓂,𝓃=0,1,\dots ,N \\ {\mathit{H}}_{f,i}=\left[\begin{array}{cccc}{h}_{00}^{f,i}& h_{01}^{f,i}& \cdots & h_{0N}^{f,i}\\ h_{10}^{f,i}& h_{11}^{f,i}& \cdots & h_{1N}^{f,i}\\ ⋮& ⋮& \ddots & ⋮\\ h_{N0}^{f,i}& h_{N1}^{f,i}& \cdots & h_{NN}^{f,i}\end{array}\right], {\mathit{H}}_{v,j}\left(x_{\mathit{𝒾}}\right)=\left[\begin{array}{cccc}{h}_{00}^{v,j}\left(x_{\mathit{𝒾}}\right)& h_{01}^{v,j}\left(x_{\mathit{𝒾}}\right)& \cdots & h_{0N}^{v,j}\left(x_{\mathit{𝒾}}\right)\\ h_{10}^{v,j}\left(x_{\mathit{𝒾}}\right)& h_{11}^{v,j}\left(x_{\mathit{𝒾}}\right)& \cdots & h_{1N}^{v,j}\left(x_{\mathit{𝒾}}\right)\\ ⋮& ⋮& \ddots & ⋮\\ h_{N0}^{v,j}\left(x_{\mathit{𝒾}}\right)& h_{N1}^{v,j}\left(x_{\mathit{𝒾}}\right)& \cdots & h_{NN}^{v,j}\left(x_{\mathit{𝒾}}\right)\end{array}\right],\mathit{𝒾}=0,1,\dots ,N \end{gathered}$$

where, respectively

$$\begin{gathered} \left[h_{rs}^{f,i}\right]=\frac{b^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i+r+s+1}-a^{\overline{n}\left({\alpha }_i\right)-{\alpha }_i+r+s+1}}{\overline{n}\left({\alpha }_i\right)-{\alpha }_i+r+s+1}, i=0,1,\dots ,m_1, r,s=0,1,\dots ,N \\ \left[h_{rs}^{v,j}\left(x_{\mathrm{𝒾}}\right)\right]=\frac{x_{\mathrm{𝒾}}{}^{\overline{n}\left({\beta }_j\right)-{\beta }_j+r+s+1}-a^{\overline{n}\left({\beta }_j\right)-{\beta }_j+r+s+1}}{\overline{n}\left({\beta }_j\right)-{\beta }_j+r+s+1}, j=0,1,\dots ,m_2, r,s=0,1,\dots ,N \\ \mathrm{𝒾}=0,1,\dots ,N \\ {\overline{\mathit{V}}}^j=\left[\begin{array}{cccc}{\mathit{V}}_t^j& 0& \cdots & 0\\ 0& {\mathit{V}}_t^j& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\mathit{V}}_t^j\end{array}\right],{\overline{\mathit{C}}}_j=\left[\begin{array}{cccc}\mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right)& 0& \cdots & 0\\ 0& \mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right)& \cdots & 0\\ ⋮& \ddots & \ddots & ⋮\\ 0& 0& \cdots & \mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right)\end{array}\right] \\ {\overline{\mathit{H}}}_j=\left[\begin{array}{cccc}{\mathit{H}}_{v,j}\left(x_0\right)& 0& \cdots & 0\\ 0& {\mathit{H}}_{v,j}\left(x_1\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\mathit{H}}_{v,j}\left(x_N\right)\end{array}\right], {\overline{\mathit{B}}}^j=\left[\begin{array}{cccc}{\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}& 0& \cdots & 0\\ 0& {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}\end{array}\right] \\ j=0,1,\dots ,m_2 \\ \overline{\mathit{D}}=\left[\begin{array}{c}{\mathit{D}}^T\\ {\mathit{D}}^T\\ ⋮\\ {\mathit{D}}^T\end{array}\right], \mathit{G}=\left[\begin{array}{c}g\left(x_0\right)\\ g\left(x_1\right)\\ ⋮\\ g\left(x_N\right)\end{array}\right], and \mathit{A}=\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_N\end{array}\right] \end{gathered}$$

When the matrices

$$\mathit{y}\left(n\right), \mathit{y}\left(n-l\right), {\mathit{p}}_n, {\mathit{p}}_l, \mathit{X}, \mathit{C}\left(\overline{n}\left({\sigma }_n\right)-{\sigma }_n\right), \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right), \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right), \mathit{C}\left(\overline{n}\left({\beta }_j\right)-{\beta }_j\right), {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)},$$

$${\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}, {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}, {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_j\right)}, {\mathit{F}}_t^i, {\mathit{V}}_t^j, {\mathit{H}}_{f,i}, {\mathit{H}}_{v,j}\left(x_{\mathrm{𝒾}}\right) \mathrm{and} {\mathit{D}}^T$$

For all $l=1,2,\dots ,n-1, i=0,1,\dots ,m_1, j=0,1,\dots ,m_2, \mathrm{𝒾}=0,1,\dots ,N$. In Equation (31) we have explained that their dimensions are similar to those of $\left(N+1\right)\times \left(N+1\right)$. Moreover, in Equation (31), these matrices $\overline{\mathit{X}}, {\overline{\mathit{V}}}^j, {\overline{\mathit{H}}}_j, {\overline{\mathit{C}}}_j, {\overline{\mathit{B}}}^j \mathrm{and} \overline{\mathit{D}}, \mathrm{for}\mathrm{all} j=0,1,\dots ,m_2$ are written in full, their measured dimensions can be observed by $\left(N+1\right)\times {\left(N+1\right)}^2,{\left(N+1\right)}^2\times {\left(N+1\right)}^2, {\left(N+1\right)}^2\times {\left(N+1\right)}^2,{\left(N+1\right)}^2\times {\left(N+1\right)}^2,{\left(N+1\right)}^2\times {\left(N+1\right)}^2$, and ${\left(N+1\right)}^2\times \left(N+1\right)$ respectively.

As a result, the fundamental matrix Equation (31) that corresponds to Equation (1) may be expressed as

$$\mathit{W}\mathit{A}=\mathit{G} or \left[\mathit{W}:\mathit{G}\right].$$

(32)

where

$$\begin{gathered} \mathit{W}=\mathit{y}\left(n\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_n\right)-{\sigma }_n\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_n\right)}{\mathit{D}}^T+{\displaystyle \sum }_{l=1}^{n-1}{\mathit{p}}_l\mathit{y}\left(n-l\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_{n-l}\right)-{\sigma }_{n-l}\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_{n-l}\right)}{\mathit{D}}^T \\ +{\mathit{p}}_n\mathit{X}{\mathit{D}}^T-{\displaystyle \sum }_{i=0}^{m_1}{\lambda }_i\mathit{X}{\mathit{F}}_t^i {\mathit{H}}_{f,i} \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T-{\displaystyle \sum }_{j=0}^{m_2}{\overline{\lambda }}_j\overline{\mathit{X}}{\overline{\mathit{V}}}^j{\overline{\mathit{H}}}_j{\overline{\mathit{C}}}_j{\overline{\mathit{B}}}^j\overline{\mathit{D}}. \end{gathered}$$

Note that, Equation (32) is a set of $\left(N+1\right)$ linear algebraic equations with unknown Bessel coefficients $\mathit{A}=\left[a_0, a_1, \dots , a_N\right].$ The matrix form (29) for the conditions, on the other hand, may be represented as

$${\mathit{U}}_k\mathit{A}=\left[{\mathit{C}}_k\right] \mathrm{or} \left[{\mathit{U}}_k:{\mathit{C}}_k\right]; k=0,1,\dots ,\mu -1, \mu =\max \left\{\lceil {\sigma }_n\rceil ,\lceil {\alpha }_{m_1}\rceil ,\lceil {\beta }_{m_2}\rceil \right\}.$$

(33)

$$\begin{gathered} {\mathit{U}}_k={\displaystyle \sum }_{\ell =0}^{\mu -1}\left[{𝒽}_{k\ell }\mathit{X}\left(a\right)+{\overline{𝒽}}_{k\ell }\mathit{X}\left(b\right)\right]{\left({\mathit{B}}^T\right)}^{\ell }{\mathit{D}}^T\mathit{A}. \\ =\left[u_{k0} u_{k1} u_{k2} \dots u_{kN}\right], k=0,1,\dots ,\mu -1 \end{gathered}$$

Hence, we may solve Equation (1) under mixed conditions (2) by substituting the rows of the matrices $\mathit{W}$ and $\mathit{G}$ for the rows of the matrices ${\mathit{U}}_k \mathrm{and} {\mathit{C}}_k,$ respectively.

$$\tilde{\mathit{W}}\mathit{A}=\tilde{\mathit{G}}.$$

The new augmented matrix (some time may be symmetry) of the preceding system is as follows if the last $\mu$-rows of the matrix (32) are replaced for simplicity:

$$\left[\tilde{\mathit{W}}:\tilde{\mathit{G}}\right]=\left[\begin{array}{ccccccc}{w}_{00}& w_{01}& w_{02}& \cdots & w_{0N}& :& g\left(x_0\right)\\ w_{10}& w_{11}& w_{12}& \cdots & w_{1N}& :& g\left(x_1\right)\\ w_{20}& w_{21}& w_{22}& \cdots & w_{2N}& :& g\left(x_2\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ w_{N-m,0}& w_{N-m,1}& w_{N-m,2}& \cdots & w_{N-m,N}& :& g\left(x_{N-m}\right)\\ u_{00}& u_{01}& u_{02}& \cdots & u_{0N}& :& c_0\\ u_{10}& u_{11}& u_{12}& \cdots & u_{1N}& :& c_1\\ u_{20}& u_{21}& u_{22}& \cdots & u_{2N}& :& c_2\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ u_{\mu -1,0}& u_{\mu -1,1}& u_{\mu -1,2}& \cdots & u_{\mu -1,N}& :& c_{\mu -1}\end{array}\right]$$

(34)

Take note that rank $\tilde{\mathit{W}}=\mathrm{rank}\left[\tilde{\mathit{W}}:\tilde{\mathit{G}}\right]=N+1.$ If it isn’t, the suggested technique fails to offer a solution; but in this case, in this situation, the number of collocation points (or, equivalently, the dimension of the matrix $\tilde{\mathrm{W}}$) can be increased to get the specific or general answer. As a result, we may write $\mathit{A}={\left(\tilde{\mathit{W}}\right)}^{-1}\tilde{\mathit{G}}$, and therefore the elements $a_0, a_1, \dots , a_N of \mathit{A}$ are uniquely determined.

Moreover, select that $N$ we define needs to be greater than $\mu$, i.e., $N>\mu =max\left\{\lceil {\sigma }_n\rceil ,\lceil {\alpha }_{m_1}\rceil , \lceil {\beta }_{m_2}\rceil \right\}.$ If it is not, the proposed strategy is thus unable to give a solution, because matrix ${\mathit{B}}^T$ becomes a zero matrix, we only get zero solution.

5. Numerical Examples

In this work, we choose several examples where the exact solution already exists to demonstrate the accuracy. They were all carried out on a computer using a Python program V3.8.8 (2021). The least square errors (L.S.E) in tables are the values of $\sum _{\mathrm{𝒾}=0}^M{\left[u\left(x_{\mathrm{𝒾}}\right)-{\tilde{u}}_N\left(x_{\mathrm{𝒾}}\right)\right]}^2, M\in \mathbb{N}$ at $M$-selected collocation points $x_{\mathrm{𝒾}}.$ and the running time is also provided in tabular form.

Example 1.

Consider the linear Fredholm–Volterra integro-differential equation of multi-higher fractional order, given by

$${}_0{}^{C}D{}_x^{1.2}u\left(x\right)+x{}_0{}^{C}D{}_x^{0.1}u\left(x\right)-xu\left(x\right)=g\left(x\right)+\underset{0}{\overset{1}{{\displaystyle \int }}}\left(e^{x}u\left(t\right)-t{}_0{}^{C}D{}_t^{0.9}u\left(t\right)\right)dt+\underset{0}{\overset{x}{{\displaystyle \int }}}\left(-\left(x-t\right)u\left(t\right)+2x{}_0{}^{C}D{}_t^{1.8}u\left(t\right)\right)dt$$

where $0\le x,t\le 1$

$$g\left(x\right)=\frac{-2}{\mathsf{\Gamma }\left(1.8\right)}{x}^{0.8}-\frac{2}{\mathsf{\Gamma }\left(2.9\right)}{x}^{2.9}+\frac{1}{\mathsf{\Gamma }\left(1.9\right)}{x}^{1.9}+\frac{4}{\mathsf{\Gamma }\left(2.2\right)}{x}^{2.2}-\frac{1}{6}{e}^x+\frac{1.1}{\mathsf{\Gamma }\left(3.1\right)}-\frac{4.2}{\mathsf{\Gamma }\left(4.1\right)}-x^2+\frac{7}{6}{x}^3-\frac{1}{12}{x}^4$$

with the boundary conditions

$$2u^{\left(1\right)}\left(0\right)+u^{\left(1\right)}\left(1\right)=1 \mathrm{and} u^{\left(1\right)}\left(1\right)=-1$$

which is the exact solution $u\left(x\right)=x\left(1-x\right).$

Let us now determine the $N$-truncated Bessel series approximate solution $u_N\left(x\right)$

$$u\left(x\right)\cong u_N\left(x\right)={\displaystyle \sum }_{𝓇=0}^N{a}_{𝓇}{\mathrm{J}}_{𝓇}\left(x\right)$$

Here, from the considered, example we have:

$$\begin{array}{cc}{\sigma }_0=0,{\sigma }_1=0.1,{\sigma }_2=1.2 \to \overline{n}\left({\sigma }_0\right)={\sigma }_0=0,\overline{n}\left({\sigma }_1\right)={\sigma }_1=1,\overline{n}\left({\sigma }_2\right)={\sigma }_2=2\hfill & \hfill \\ {\alpha }_0=0,{\alpha }_1=0.9 \to \overline{n}\left({\alpha }_0\right)={\alpha }_0=0,\overline{n}\left({\alpha }_1\right)={\alpha }_1=1\hfill & \hfill \\ {\beta }_0=0,{\beta }_1=1.8 \to \overline{n}\left({\beta }_0\right)={\beta }_0=0,\overline{n}\left({\beta }_1\right)={\beta }_1=2\hfill & \hfill \\ \mu =max\left\{\lceil 1.2\rceil ,\lceil 0.9\rceil , \lceil 1.8\rceil \right\}=2\hfill & \hfill \\ {\mathit{p}}_1\left(x\right)=x,{\mathit{p}}_2\left(x\right)=-x,F_0\left(x,t\right)=e^x,F_1\left(x,t\right)=t,V_0\left(x,t\right)=x-t,V_1\left(x,t\right)=x,\hfill & \hfill \\ {\lambda }_0=1,{\lambda }_1=-1,{\overline{\lambda }}_0=-1,{\overline{\lambda }}_1=2\hfill & \hfill \end{array}$$

Hence $\mu =2$ so take, the collocation point sets are $\left\{x_0=0, x_1=\frac{1}{3}, x_2=\frac{2}{3}, x_3=1\right\},$ and the fundamental matrix equation of the given (LF-VFIDEs) is derived from Equation (31), written as

$$\begin{gathered} [\mathit{y}\left(2\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_2\right)-{\sigma }_2\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_2\right)}{\mathit{D}}^T+{\mathit{p}}_1\mathit{y}\left(1\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_1\right)-{\sigma }_1\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_1\right)}{\mathit{D}}^T+{\mathit{p}}_2\mathit{X}{\mathit{D}}^T \\ -{\displaystyle \sum }_{i=0}^1{\lambda }_i\mathit{X}{\mathit{F}}_t^i {\mathit{H}}_{f,i} \mathit{C}\left(\overline{n}\left({\alpha }_i\right)-{\alpha }_i\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_i\right)}{\mathit{D}}^T-{\displaystyle \sum }_{j=0}^1{\overline{\lambda }}_j\overline{\mathit{X}}{\overline{\mathit{V}}}^j{\overline{\mathit{H}}}_j{\overline{\mathit{C}}}_j{\overline{\mathit{B}}}^j\overline{\mathit{D}} ]\mathit{A}=\mathit{G} \end{gathered}$$

where

$$\begin{gathered} \mathit{X}\left(x\right)=\left[1 x x^2 x^3\right] ,\mathit{X}=\left[\begin{array}{c}\mathit{X}\left(0\right)\\ \mathit{X}\left(1/3\right)\\ \mathit{X}\left(2/3\right)\\ \mathit{X}\left(1\right)\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 1/3& 1/9& 1/27\\ 1& 2/3& 4/9& 8/27\\ 1& 1& 1& 1\end{array}\right], \\ {\mathit{p}}_1=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1/3& 0& 0\\ 0& 0& 2/3& 0\\ 0& 0& 0& 1\end{array}\right],{\mathit{p}}_2=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& -\frac{1}{3}& 0& 0\\ 0& 0& -\frac{2}{3}& 0\\ 0& 0& 0& -1\end{array}\right], \\ \mathit{y}\left(2\right)=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& {\left(1/3\right)}^{0.8}& 0& 0\\ 0& 0& {\left(2/3\right)}^{0.8}& 0\\ 0& 0& 0& 1\end{array}\right],\mathit{y}\left(1\right)=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& {\left(1/3\right)}^{0.9}& 0& 0\\ 0& 0& {\left(2/3\right)}^{0.9}& 0\\ 0& 0& 0& 1\end{array}\right] \\ \mathit{C}\left(\overline{n}\left({\sigma }_2\right)-{\sigma }_2\right)=\left[\begin{array}{cccc}\frac{1}{\mathsf{\Gamma }\left(1.8\right)}& 0& 0& 0\\ 0& \frac{1}{\mathsf{\Gamma }\left(2.8\right)}& 0& 0\\ 0& 0& \frac{2}{\mathsf{\Gamma }\left(3.8\right)}& 0\\ 0& 0& 0& \frac{6}{\mathsf{\Gamma }\left(4.8\right)}\end{array}\right],\mathit{C}\left(\overline{n}\left({\sigma }_1\right)-{\sigma }_1\right)=\left[\begin{array}{cccc}\frac{1}{\mathsf{\Gamma }\left(1.9\right)}& 0& 0& 0\\ 0& \frac{1}{\mathsf{\Gamma }\left(2.9\right)}& 0& 0\\ 0& 0& \frac{2}{\mathsf{\Gamma }\left(3.9\right)}& 0\\ 0& 0& 0& \frac{6}{\mathsf{\Gamma }\left(4.9\right)}\end{array}\right], \\ \mathit{C}\left(\overline{n}\left({\alpha }_1\right)-{\alpha }_1\right)=\left[\begin{array}{cccc}\frac{1}{\mathsf{\Gamma }\left(1.1\right)}& 0& 0& 0\\ 0& \frac{1}{\mathsf{\Gamma }\left(2.1\right)}& 0& 0\\ 0& 0& \frac{2}{\mathsf{\Gamma }\left(3.1\right)}& 0\\ 0& 0& 0& \frac{6}{\mathsf{\Gamma }\left(4.1\right)}\end{array}\right],\mathit{C}\left(\overline{n}\left({\beta }_1\right)-{\beta }_1\right)=\left[\begin{array}{cccc}\frac{1}{\mathsf{\Gamma }\left(1.2\right)}& 0& 0& 0\\ 0& \frac{1}{\mathsf{\Gamma }\left(2.2\right)}& 0& 0\\ 0& 0& \frac{2}{\mathsf{\Gamma }\left(3.2\right)}& 0\\ 0& 0& 0& \frac{6}{\mathsf{\Gamma }\left(4.2\right)}\end{array}\right] \\ {\mathit{B}}^T={\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_1\right)}={\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_1\right)}=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 2& 0\\ 0& 0& 0& 3\\ 0& 0& 0& 0\end{array}\right], {\left({\mathit{B}}^T\right)}^2={\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_2\right)}={\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_1\right)}=\left[\begin{array}{cccc}0& 0& 2& 0\\ 0& 0& 0& 6\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right], \\ {\mathit{F}}_t^0=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& 0& 0\\ \frac{1}{2}& 0& 0& 0\\ \frac{1}{6}& 0& 0& 0\end{array}\right], {\mathit{F}}_t^1=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],{\mathit{V}}_t^0=\left[\begin{array}{cccc}0& -1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right], {\mathit{V}}_t^1=\left[\begin{array}{cccc}0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right], \\ {\mathit{H}}_{f,0}=\left[\begin{array}{cccc}\frac{1}{1}& \frac{1}{2}& \frac{1}{3}& \frac{1}{4}\\ \frac{1}{2}& \frac{1}{3}& \frac{1}{4}& \frac{1}{5}\\ \frac{1}{3}& \frac{1}{4}& \frac{1}{5}& \frac{1}{6}\\ \frac{1}{4}& \frac{1}{5}& \frac{1}{6}& \frac{1}{7}\end{array}\right], {\mathit{H}}_{f,1}=\left[\begin{array}{cccc}\frac{1}{1.1}& \frac{1}{2.1}& \frac{1}{3.1}& \frac{1}{4.1}\\ \frac{1}{2.1}& \frac{1}{3.1}& \frac{1}{4.1}& \frac{1}{5.1}\\ \frac{1}{3.1}& \frac{1}{4.1}& \frac{1}{5.1}& \frac{1}{6.1}\\ \frac{1}{4.1}& \frac{1}{5.1}& \frac{1}{6.1}& \frac{1}{7.1}\end{array}\right], {\mathit{D}}^T=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \frac{1}{2}& 0& 0\\ \frac{-1}{4}& 0& \frac{1}{8}& 0\\ 0& \frac{-1}{16}& 0& \frac{1}{48}\end{array}\right], \\ \overline{\mathit{X}}=\left[\begin{array}{cccc}\mathit{X}\left(0\right)& 0& 0& 0\\ 0& \mathit{X}\left(1/3\right)& 0& 0\\ 0& 0& \mathit{X}\left(2/3\right)& 0\\ 0& 0& 0& \mathit{X}\left(1\right)\end{array}\right],{\overline{\mathit{V}}}^0=\left[\begin{array}{cccc}{\mathit{V}}_t^0& 0& 0& 0\\ 0& {\mathit{V}}_t^0& 0& 0\\ 0& 0& {\mathit{V}}_t^0& 0\\ 0& 0& 0& {\mathit{V}}_t^0\end{array}\right], {\overline{\mathit{V}}}^1=\left[\begin{array}{cccc}{\mathit{V}}_t^1& 0& 0& 0\\ 0& {\mathit{V}}_t^1& 0& 0\\ 0& 0& {\mathit{V}}_t^1& 0\\ 0& 0& 0& {\mathit{V}}_t^1\end{array}\right] \\ {\overline{\mathit{H}}}_0=\left[\begin{array}{cccc}{\mathit{H}}_{v,0}\left(0\right)& 0& 0& 0\\ 0& {\mathit{H}}_{v,0}\left(1/3\right)& 0& 0\\ 0& 0& {\mathit{H}}_{v,0}\left(2/3\right)& 0\\ 0& 0& 0& {\mathit{H}}_{v,0}\left(1\right)\end{array}\right],{\overline{\mathit{H}}}_1=\left[\begin{array}{cccc}{\mathit{H}}_{v,1}\left(0\right)& 0& 0& 0\\ 0& {\mathit{H}}_{v,1}\left(1/3\right)& 0& 0\\ 0& 0& {\mathit{H}}_{v,1}\left(2/3\right)& 0\\ 0& 0& 0& {\mathit{H}}_{v,1}\left(1\right)\end{array}\right], \\ {\mathit{H}}_{v,0}\left(0\right)=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],{\mathit{H}}_{v,0}\left(1/3\right)=\left[\begin{array}{cccc}\frac{1}{3}& \frac{{\left(\frac{1}{3}\right)}^2}{2}& \frac{{\left(\frac{1}{3}\right)}^3}{3}& \frac{{\left(\frac{1}{3}\right)}^4}{4}\\ \frac{{\left(\frac{1}{3}\right)}^2}{2}& \frac{{\left(\frac{1}{3}\right)}^3}{3}& \frac{{\left(\frac{1}{3}\right)}^4}{4}& \frac{{\left(\frac{1}{3}\right)}^5}{5}\\ \frac{{\left(\frac{1}{3}\right)}^3}{3}& \frac{{\left(\frac{1}{3}\right)}^4}{4}& \frac{{\left(\frac{1}{3}\right)}^5}{5}& \frac{{\left(\frac{1}{3}\right)}^6}{6}\\ \frac{{\left(\frac{1}{3}\right)}^4}{4}& \frac{{\left(\frac{1}{3}\right)}^5}{5}& \frac{{\left(\frac{1}{3}\right)}^6}{6}& \frac{{\left(\frac{1}{3}\right)}^7}{7}\end{array}\right], {\mathit{H}}_{v,0}\left(2/3\right)=\left[\begin{array}{cccc}\frac{2}{3}& \frac{{\left(\frac{2}{3}\right)}^2}{2}& \frac{{\left(\frac{2}{3}\right)}^3}{3}& \frac{{\left(\frac{2}{3}\right)}^4}{4}\\ \frac{{\left(\frac{2}{3}\right)}^2}{2}& \frac{{\left(\frac{2}{3}\right)}^3}{3}& \frac{{\left(\frac{2}{3}\right)}^4}{4}& \frac{{\left(\frac{2}{3}\right)}^5}{5}\\ \frac{{\left(\frac{2}{3}\right)}^3}{3}& \frac{{\left(\frac{2}{3}\right)}^4}{4}& \frac{{\left(\frac{2}{3}\right)}^5}{5}& \frac{{\left(\frac{2}{3}\right)}^6}{6}\\ \frac{{\left(\frac{2}{3}\right)}^4}{4}& \frac{{\left(\frac{2}{3}\right)}^5}{5}& \frac{{\left(\frac{2}{3}\right)}^6}{6}& \frac{{\left(\frac{2}{3}\right)}^7}{7}\end{array}\right], \\ {\mathit{H}}_{v,0}\left(1\right)=\left[\begin{array}{cccc}\frac{1}{1}& \frac{1}{2}& \frac{1}{3}& \frac{1}{4}\\ \frac{1}{2}& \frac{1}{3}& \frac{1}{4}& \frac{1}{5}\\ \frac{1}{3}& \frac{1}{4}& \frac{1}{5}& \frac{1}{6}\\ \frac{1}{4}& \frac{1}{5}& \frac{1}{6}& \frac{1}{7}\end{array}\right],{\mathit{H}}_{v,1}\left(0\right)=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],{\mathit{H}}_{v,1}\left(1/3\right)=\left[\begin{array}{cccc}\frac{{\left(\frac{1}{3}\right)}^{1.2}}{1.2}& \frac{{\left(\frac{1}{3}\right)}^{2.2}}{2.2}& \frac{{\left(\frac{1}{3}\right)}^{3.2}}{3.2}& \frac{{\left(\frac{1}{3}\right)}^{4.2}}{4.2}\\ \frac{{\left(\frac{1}{3}\right)}^{2.2}}{2.2}& \frac{{\left(\frac{1}{3}\right)}^{3.2}}{3.2}& \frac{{\left(\frac{1}{3}\right)}^{4.2}}{4.2}& \frac{{\left(\frac{1}{3}\right)}^{5.2}}{5.2}\\ \frac{{\left(\frac{1}{3}\right)}^{3.2}}{3.2}& \frac{{\left(\frac{1}{3}\right)}^{4.2}}{4.2}& \frac{{\left(\frac{1}{3}\right)}^{5.2}}{5.2}& \frac{{\left(\frac{1}{3}\right)}^{6.2}}{6.2}\\ \frac{{\left(\frac{1}{3}\right)}^{4.2}}{4.2}& \frac{{\left(\frac{1}{3}\right)}^{5.2}}{5.2}& \frac{{\left(\frac{1}{3}\right)}^{6.2}}{6.2}& \frac{{\left(\frac{1}{3}\right)}^{7.2}}{7.2}\end{array}\right], \\ {\mathit{H}}_{v,1}\left(2/3\right)=\left[\begin{array}{cccc}\frac{{\left(\frac{2}{3}\right)}^{1.2}}{1.2}& \frac{{\left(\frac{2}{3}\right)}^{2.2}}{2.2}& \frac{{\left(\frac{2}{3}\right)}^{3.2}}{3.2}& \frac{{\left(\frac{2}{3}\right)}^{4.2}}{4.2}\\ \frac{{\left(\frac{2}{3}\right)}^{2.2}}{2.2}& \frac{{\left(\frac{2}{3}\right)}^{3.2}}{3.2}& \frac{{\left(\frac{2}{3}\right)}^{4.2}}{4.2}& \frac{{\left(\frac{2}{3}\right)}^{5.2}}{5.2}\\ \frac{{\left(\frac{2}{3}\right)}^{3.2}}{3.2}& \frac{{\left(\frac{2}{3}\right)}^{4.2}}{4.2}& \frac{{\left(\frac{2}{3}\right)}^{5.2}}{5.2}& \frac{{\left(\frac{2}{3}\right)}^{6.2}}{6.2}\\ \frac{{\left(\frac{2}{3}\right)}^{4.2}}{4.2}& \frac{{\left(\frac{2}{3}\right)}^{5.2}}{5.2}& \frac{{\left(\frac{2}{3}\right)}^{6.2}}{6.2}& \frac{{\left(\frac{2}{3}\right)}^{7.2}}{7.2}\end{array}\right],{\mathit{H}}_{v,1}\left(1\right)=\left[\begin{array}{cccc}\frac{1}{1.2}& \frac{1}{2.2}& \frac{1}{3.2}& \frac{1}{4.2}\\ \frac{1}{2.2}& \frac{1}{3.2}& \frac{1}{4.2}& \frac{1}{5.2}\\ \frac{1}{3.2}& \frac{1}{4.2}& \frac{1}{5.2}& \frac{1}{6.2}\\ \frac{1}{4.2}& \frac{1}{5.2}& \frac{1}{6.2}& \frac{1}{7.2}\end{array}\right], \\ {\overline{\mathit{C}}}_1=\left[\begin{array}{cccc}\mathit{C}\left(\overline{n}\left({\beta }_1\right)-{\beta }_1\right)& 0& 0& 0\\ 0& \mathit{C}\left(\overline{n}\left({\beta }_1\right)-{\beta }_1\right)& 0& 0\\ 0& 0& \mathit{C}\left(\overline{n}\left({\beta }_1\right)-{\beta }_1\right)& 0\\ 0& 0& 0& \mathit{C}\left(\overline{n}\left({\beta }_1\right)-{\beta }_1\right)\end{array}\right], \\ {\overline{\mathit{B}}}^1=\left[\begin{array}{cccc}{\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_1\right)}& 0& 0& 0\\ 0& {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_1\right)}& 0& 0\\ 0& 0& {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_1\right)}& 0\\ 0& 0& 0& {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\beta }_1\right)}\end{array}\right], \mathit{G}=\left[\begin{array}{c}-0.28262788 \\ -0.90165393\\ -0.47693449\\ 0.94267344\end{array}\right] and \overline{\mathit{D}}=\left[\begin{array}{c}{\mathit{D}}^T\\ {\mathit{D}}^T\\ {\mathit{D}}^T\\ {\mathit{D}}^T\end{array}\right] \end{gathered}$$

putting all above matrices in matrix Equation (32) and calculating it, this fundamental matrix equation’s augmented matrix is:

$$\left[\mathit{W}:\mathit{G}\right]=\left[\begin{array}{cccccc}-1.07079233& -0.02572358& 0.03539616& 0.00866477& :& -0.28262788\\ -1.85498405& -0.12829141& 0.09107227& 0.01393314& :& -0.90165393\\ -2.40595015& -0.22844205& 0.01161705& 0.01241051& :& -0.47693449\\ -2.7722549& -0.23879802& -0.19720588& -0.02479584& :& 0.94267344\end{array}\right]$$

For our consider example, the boundary conditions from Equation (33) have the following matrix forms:

$${\mathit{U}}_k\mathit{A}=\left[{\mathit{C}}_k\right] or \left[{\mathit{U}}_k:{\mathit{C}}_k\right]; k=0,1$$

or clearly

$$\begin{array}{cc}\left[{\mathit{U}}_0:{\mathit{C}}_0\right]=\left[\begin{array}{cccccc}-0.5& 1.3125& 0.25& 0.0625& :& 1\end{array}\right]\hfill & \hfill \\ \left[{\mathit{U}}_1:{\mathit{C}}_1\right]=\left[\begin{array}{cccccc}-0.5& 0.3125& 0.25& 0.0625& :& -1\end{array}\right] \hfill & \hfill \end{array}$$

The new augmented matrix depending on conditions is constructed as follows from the system (34):

$$\left[\tilde{\mathit{W}}:\tilde{\mathit{G}}\right]=\left[\begin{array}{cccccc}-1.07079233& -0.02572358& 0.03539616& 0.00866477& :& -0.28262788\\ -1.85498405& -0.12829141& 0.09107227& 0.01393314& :& -0.90165393\\ -0.5& 1.3125& 0.25& 0.0625& :& 1\\ -0.5& 0.3125& 0.25& 0.0625& :& -1\end{array}\right]$$

The Bessel coefficient matrix $\mathit{A}$ is obtained by solving this system.

$$\mathit{A}={\left[-1.98357176\times {10}^{-06} 2.00000000 -8.00269502 6.01076421\right]}^T$$

hence, for $N=3$ the approximate solution of the problem is formed as

$$u_3\left(x\right)=0.0002242544x^3-1.0003363816x^2+1.0x-0.0000019836$$

for $N=6$ and $N=10$, similarly as steps above and running the general python program which are written for this purpose we obtain the approximate solution of the problem, respectively.

$$\begin{gathered} u_6\left(x\right)=-0.0003537963x^6+0.0007932088x^5-0.0006209581x^4+0.000229553x^3 \\ -1.0000240462x^2+1.0x+0.0000090551 \end{gathered}$$

and

$$\begin{gathered} u_{10}\left(x\right)=0.0000000661x^{10}-0.0000002889x^9+0.00000052947x^8-0.00000053514x^7 \\ \begin{array}{cc}\hfill & +0.000000328x^6-0.00000012512x^5+0.0000000293x^4\hfill \\ \hfill & -0.00000000370x^3-0.999999999656939x^2+1.0x^{ }+0.00000000015\hfill \end{array} \end{gathered}$$

In

Table 1. Compares the exact solution $u\left(x\right)$ with the approximate solution $u_N\left(x\right)$ of Example 1.

${\mathit{x}}_{\mathbf{𝒾}}$	Exact Solution Example 1.	$\mathit{N}$$\text{-}\mathbf{Approximate}\mathbf{Solution}{\mathit{u}}_{\mathit{N}}\left(\mathit{x}\right)$
		N = 3	$\mathit{N}=\mathbf{6}$	$\mathit{N}=\mathbf{10}$
$0.0$	$0.00$	$-1.9835710\times {10}^{-06}$	$9.05511829\times {10}^{-06}$	$1.53060853\times {10}^{-10}$
$0.1$	$0.09$	$8.99948769\times {10}^{-02}$	9.00089897 × 10−02	$9.00000002\times {10}^{-02}$
$0.2$	$0.16$	$1.59986355\times {10}^{-01}$	$1.60009167\times {10}^{-01}$	$1.60000000\times {10}^{-01}$
$0.3$	$0.21$	2.09973797 × 10−01	$2.10009729\times {10}^{-01}$	$2.10000000\times {10}^{-01}$
$0.4$	$0.24$	$2.39958548\times {10}^{-01}$	$2.40010676\times {10}^{-01}$	$2.40000000\times {10}^{-01}$
$0.5$	$0.25$	$2.49941953\times {10}^{-01}$	$2.50012188\times {10}^{-01}$	$2.50000000\times {10}^{-01}$
$0.6$	$0.24$	$2.39925358\times {10}^{-01}$	$2.40014679\times {10}^{-01}$	$2.40000000\times {10}^{-01}$
$0.7$	$0.21$	$2.09910109\times {10}^{-01}$	$2.10018608\times {10}^{-01}$	$2.10000000\times {10}^{-01}$
$0.8$	$0.16$	$1.59897550\times {10}^{-01}$	$1.60024025\times {10}^{-01}$	$1.60000000\times {10}^{-01}$
$0.9$	$0.09$	$8.98890288\times {10}^{-02}$	$9.00298712\times {10}^{-02}$	$9.00000005\times {10}^{-02}$
$1.0$	$0.00$	$-1.141107\times {10}^{-04}$	$3.30162295\times {10}^{-05}$	$5.74892312\times {10}^{-10}$
L.S.E.		$5.5474382\times {10}^{-08}$	$3.72530728\times {10}^{-09}$	$1.13502895\times {10}^{-18}$
Running Time/Sec		$0.47589683$	$0.227055311$	$0.57032561$

, Comparison the exact solution $u\left(x\right)$ with the approximate solution $u_N\left(x\right)$ of Example 1 for $N=3,6\mathrm{and}10$, respectively, in terms of least square error and running time.

Example 2.

Let us now consider the LF-VIFDEs on the closed bounded interval $\left[0,1\right]$ given by

$$\begin{array}{cc}\hfill {}_0{}^{C}D{}_x^{1.3}u\left(x\right)+\frac{x}{2}{}_0{}^{C}D{}_x^{0.8}u\left(x\right)& +\sqrt{x}{}_0{}^{C}D{}_x^{0.5}u\left(x\right)+\left(x^2+1\right)u\left(x\right)\hfill \\ \hfill =g\left(x\right)& +\underset{0}{\overset{1}{{\displaystyle \int }}}\left(\left(\sin \left(x\right)-t\right){}_0{}^{C}D{}_t^{0.9}u\left(t\right)+2\left(e^x-t\right){}_0{}^{C}D{}_t^{1.1}u\left(t\right)\right)dt+3\underset{0}{\overset{x}{{\displaystyle \int }}}t\sinh \left(x\right){}_0{}^{C}D{}_t^{1.9}u\left(t\right)dt\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill g\left(x\right)=\frac{2}{\mathsf{\Gamma }\left(1.7\right)}{x}^{0.7}& +\frac{1}{\mathsf{\Gamma }\left(2.2\right)}{x}^{2.2}+\frac{1}{\mathsf{\Gamma }\left(1.2\right)}{x}^{1.2}+\frac{2}{\mathsf{\Gamma }\left(2.5\right)}{x}^{1.5}\sqrt{x}+\frac{2}{\mathsf{\Gamma }\left(1.5\right)}{x}^{ }+\left(x^2+1\right){\left(x+1\right)}^2\hfill \\ & -\left(\frac{2}{\mathsf{\Gamma }\left(3.1\right)}+\frac{2}{\mathsf{\Gamma }\left(2.1\right)}\right)\sin \left(x\right)-\frac{4}{\mathsf{\Gamma }\left(2.9\right)}{e}^x-\frac{6.6}{\mathsf{\Gamma }\left(3.1\right)}{x}^{2.1}\sinh \left(x\right)+\left(\frac{4.2}{\mathsf{\Gamma }\left(4.1\right)}+\frac{2.2}{\mathsf{\Gamma }\left(3.1\right)}+\frac{7.6}{\mathsf{\Gamma }\left(3.9\right)}\right)\hfill \end{array}$$

with the boundary conditions:

$$u^{ }\left(0\right)+u^{\left(1\right)}\left(1\right)=5 \mathrm{and} u^{ }\left(1\right)+u^{\left(1\right)}\left(0\right)=6$$

The exact solution is $u\left(x\right)={\left(x+1\right)}^2.$

Let us now calculate the coefficients $a_{𝓇}, 𝓇=\overline{0:N}$ of approximate solution with the aid of the truncated Bessel series:

$$u\left(x\right)\cong u_N\left(x\right)={\displaystyle \sum }_{𝓇=0}^N{a}_{𝓇}{\mathrm{J}}_{𝓇}\left(x\right)$$

Here, from the considered example we have:

$$\begin{array}{cc}{\sigma }_1=0.5,{\sigma }_2=0.8,,{\sigma }_3=1.3\to \overline{n}\left({\sigma }_1\right)={\sigma }_1=1,\overline{n}\left({\sigma }_2\right)={\sigma }_2=1,\overline{n}\left({\sigma }_3\right)={\sigma }_3=2\hfill & \hfill \\ {\alpha }_1=0.9,{\alpha }_2=1.1 \to \overline{n}\left({\alpha }_1\right)={\alpha }_1=1,\overline{n}\left({\alpha }_2\right)={\alpha }_2=2\hfill & \hfill \\ {\beta }_1=1.9 \to \overline{n}\left({\beta }_1\right)={\beta }_1=2,\mu =\max \left\{\lceil 1.3\rceil ,\lceil 1.1\rceil , \lceil 1.9\rceil \right\}=2 \mathrm{and} {\lambda }_1=1,{\lambda }_2=2,{\overline{\lambda }}_1=3\hfill & \hfill \\ {\mathit{p}}_1\left(x\right)=\frac{x}{2}, {\mathit{p}}_2\left(x\right)=\sqrt{x},{\mathit{p}}_3\left(x\right)=\left(x^2+1\right),F_1\left(x,t\right)=\sin \left(x\right)-t,F_2\left(x,t\right)=\left(e^x-t\right), V_1\left(x,t\right)=t\sinh \left(x\right)\hfill \end{array}$$

Hence $\mu =2$ so take $N=3$, the set of collocation points, and the fundamental matrix equation of the given (LF-VFIDEs) is derived from Equation (31), written as

$$\begin{gathered} [\mathit{y}\left(3\right)\mathit{X} \mathit{C}(\overline{n}\left({\sigma }_3\right)-{\sigma }_3) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_3\right)}{\mathit{D}}^T+{\mathit{p}}_1\mathit{y}\left(2\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_2\right)-{\sigma }_2\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_2\right)}{\mathit{D}}^T \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}& +{\mathit{p}}_2\mathit{y}\left(1\right)\mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_1\right)-{\sigma }_1\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_1\right)}{\mathit{D}}^T+{\mathit{p}}_3\mathit{X}{\mathit{D}}^T\hfill \\ & -{\lambda }_1\mathit{X}{\mathit{F}}_t^1 {\mathit{H}}_{f,1} \mathit{C}\left(\overline{n}\left({\alpha }_1\right)-{\alpha }_1\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_1\right)}{\mathit{D}}^T\hfill \\ & -{\lambda }_2\mathit{X}{\mathit{F}}_t^2 {\mathit{H}}_{f,2} \mathit{C}\left(\overline{n}\left({\alpha }_2\right)-{\alpha }_2\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_2\right)}{\mathit{D}}^T-{\overline{\lambda }}_1\overline{\mathit{X}}{\overline{\mathit{V}}}^1{\overline{\mathit{H}}}_1{\overline{\mathit{C}}}_1{\overline{\mathit{B}}}^1\overline{\mathit{D}} ]\mathit{A}=\mathit{G}\hfill \end{array} \end{gathered}$$

After inputting each of the parameters above by running the general python program, which are written for this purpose for $N=3$, we obtain the approximate solution of the problem,

$$u_3\left(x\right)=0.0029061023x^3+0.99568824199x^2+2.0015004466x^{ }+0.9984047624$$

Similarly, the approximate solution of the problem for $N=6 \mathrm{and} 11,$ respectively, we obtain

$$\begin{gathered} u_6\left(x\right)=4.353138854\times {10}^{-7}{x}^6-4.54065912\times {10}^{-5}{x}^5+5.511176789\times {10}^{-5}{x}^4 \\ \begin{array}{cc}\hfill & -1.135105524\times {10}^{-5}{x}^3+1.0000204x^2+1.999983583x\hfill \\ \hfill & +1.000013642\hfill \end{array} \end{gathered}$$

and

$$\begin{gathered} u_{11}\left(x\right)=3.2299188319\times {10}^{-7}{x}^{11}-1.5377530537\times {10}^{-6}{x}^{10} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}\hfill & +3.1687057159\times {10}^{-6}{x}^9-3.716458436\times {10}^{-6}{x}^8 \hfill \\ \hfill & +2.73904171\times {10}^{-6}{x}^7-1.319310342\times {10}^{-6}{x}^6\hfill \\ \hfill & +4.181528996\times {10}^{-7}{x}^5-8.549654914\times {10}^{-8}{x}^4 \hfill \\ \hfill & +1.082189316\times {10}^{-8}{x}^3+0.9999999996x^2+1.9999999996x^{ } \hfill \\ \hfill & +1.0000000003.\hfill \end{array} \end{gathered}$$

In

Table 2. Compares the exact solution $u\left(x\right)$ with the approximate solution $u_N\left(x\right)$ of Example 2.

${\mathit{x}}_{\mathbf{𝒾}}$	Exact Solution Example 2.	$\mathit{N}$$\text{-}\mathbf{Approximate}\mathbf{Solutions}{\mathit{u}}_{\mathit{N}}\left(\mathit{x}\right)$
		$\mathit{N}=3$	$\mathit{N}=6$	$\mathit{N}=11$
$0.0$	$1.00$	$0.99840476$	$1.00001364$	$1.00$
$0.1$	$1.21$	$1.2085146$	$1.2100122$	$1.21$
$0.2$	$1.44$	$1.43855563$	$1.44001116$	$1.44$
$0.3$	$1.69$	$1.6885453$	$1.69001058$	1.69
$0.4$	$1.96$	$1.95850105$	$1.96001056$	$1.96$
$0.5$	$2.25$	$2.24844031$	$2.25001115$	$2.25$
$0.6$	2.56	$2.55838052$	$2.56001232$	$2.56$
0.7	$2.89$	$2.88833911$	$2.89001391$	$2.89$
$0.8$	$3.24$	$3.23833352$	$3.24001556$	$3.24$
$0.9$	$3.61$	$3.60838119$	3.6100167	$3.61$
$1.0$	$4.00$	$3.99849955$	$4.00001642$	$4.00$
L.S.E.		$2.66633874\times {10}^{-05}$	$1.94276005\times {10}^{-09}$	$7.33991398\times {10}^{-19}$
Running Time/Sec			$0.2821376323$	$0.6396300792$

comparison in terms of least square error and running time the exact solution $u\left(x\right)$ with the approximate solution $u_N\left(x\right)$ of example 2 for $N=3,6\mathrm{and}11$, respectively.

Figure 1 and Figure 2 illustrate a comparison between the exact solution and approximate solution of LF-VIFDEs of Examples 1 and 2, respectively. To show the result of the proposed method to an exact solution, we present Table 1 and Table 2, respectively. Each of the plots is drawn with our Python program version 3.8.8 (2021).

Example 3.

Let us consider the linear Fredholm–Volterra fractional integro-differential equation on the closed bounded interval [0,1]:

$${}_0{}^{C}D{}_x^{0.73}u\left(x\right)=g\left(x\right)+\underset{0}{\overset{1}{{\displaystyle \int }}}\left(t \sin (\frac{x}{2}){}_0{}^{C}D{}_t^{0.64}u(t)+2(x-\frac{t}{2}){}_0{}^{C}D{}_t^{1.64}u\left(t\right)\right)dt+\underset{0}{\overset{x}{{\displaystyle \int }}}\left(-e^x{}_0{}^{C}D{}_t^{0.5}u\left(t\right)+\left(-2\right)t{e}^x{}_0{}^{C}D{}_t^{1.5}u\left(t\right)\right)dt$$

where

$$g\left(x\right)={\displaystyle \sum }_{i=0}^{\infty }\left[\frac{2x}{\mathsf{\Gamma }\left(i+2.36\right)}-\frac{x^{i+0.27}}{\mathsf{\Gamma }\left(\mathrm{i}+1.27\right)}+\frac{\sin \left(\frac{x}{2}\right)\left(i+1.36\right)}{\mathsf{\Gamma }\left(i+3.36\right)}-\frac{\left(i+1.36\right)}{\mathsf{\Gamma }\left(i+3.36\right)}-\frac{e^x{x}^{i+1.5}}{\mathsf{\Gamma }\left(\mathrm{i}+2.5\right)}-\frac{2\left(i+1.5\right)e^x{x}^{i+2.5}}{\mathsf{\Gamma }\left(\mathrm{i}+3.5\right)}\right]$$

with the boundary conditions

$$u\left(0\right)+u^{\left(1\right)}\left(0\right)=-1 \mathrm{and} u\left(1\right)-u^{\left(1\right)}\left(1\right)=1$$

which is the exact solution $u\left(x\right)=1-e^x.$

Let us now calculate the coefficients $a_{𝓇}, 𝓇=\overline{0:N}$ of approximate solution with the aid of the truncated Bessel series:

$$u\left(x\right)\cong u_N\left(x\right)={\displaystyle \sum }_{𝓇=0}^N{a}_{𝓇}{\mathrm{J}}_{𝓇}\left(x\right)$$

Here, from consider example we have:

$$\begin{array}{cc}{\sigma }_1=0.73\to \overline{n}\left({\sigma }_1\right)={\sigma }_1=1\hfill & \hfill \\ {\alpha }_1=0.64 \to \overline{n}\left({\alpha }_1\right)={\alpha }_1=1, {\alpha }_2=1.64 \to \overline{n}\left({\alpha }_2\right)={\alpha }_2=2\hfill & \hfill \\ {\beta }_1=0.5 \to \overline{n}\left({\beta }_1\right)={\beta }_1=1, {\beta }_2=1.5 \to \overline{n}\left({\beta }_1\right)={\beta }_1=2\hfill & \hfill \\ \mu =\max \left\{\lceil {\sigma }_1\rceil ,\lceil {\alpha }_2\rceil ,\lceil {\beta }_2\rceil \right\}=\max \left\{\lceil 0.73\rceil ,\lceil 1.64\rceil ,\lceil 1.5\rceil \right\}=2\hfill & \hfill \\ {\mathit{p}}_1\left(x\right)=0, F_1\left(x,t\right)=t\sin \left(\frac{x}{2}\right), F_2\left(x,t\right)=\left(x-\frac{t}{2}\right), V_1\left(x,t\right)=e^x, V_2\left(x,t\right)=t{e}^x\mathrm{and} {\lambda }_1=1,{\lambda }_2=\hfill & \hfill \\ 1, {\overline{\lambda }}_1=-1, {\overline{\lambda }}_2=-2.\hfill \end{array}$$

Suppose that, we take $\overline{N}$ terms from the homogeneous part $g\left(x\right)$:

$$g\left(x\right)={\displaystyle \sum }_{i=0}^{\overline{N}}\left[\frac{2x}{\mathsf{\Gamma }\left(i+2.36\right)}-\frac{x^{i+0.27}}{\mathsf{\Gamma }\left(\mathrm{i}+1.27\right)}+\frac{\sin \left(\frac{x}{2}\right)\left(i+1.36\right)}{\mathsf{\Gamma }\left(i+3.36\right)}-\frac{\left(i+1.36\right)}{\mathsf{\Gamma }\left(i+3.36\right)}-\frac{e^x{x}^{i+1.5}}{\mathsf{\Gamma }\left(\mathrm{i}+2.5\right)}-\frac{2\left(i+1.5\right)e^x{x}^{i+2.5}}{\mathsf{\Gamma }\left(\mathrm{i}+3.5\right)}\right]$$

Hence $\mu =2$, the fundamental matrix equation of the given (LF-VFIDEs) is derived from Equation (31), written as

$$\begin{gathered} [\mathit{y}\left(1\right)\mathit{X} \mathit{C}(\overline{n}\left({\sigma }_1\right)-{\sigma }_1) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_1\right)}{\mathit{D}}^T-{\lambda }_1\mathit{X}{\mathit{F}}_t^1 {\mathit{H}}_{f,1} \mathit{C}\left(\overline{n}\left({\alpha }_1\right)-{\alpha }_1\right){\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_1\right)}{\mathit{D}}^T \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}& -{\lambda }_2\mathit{X}{\mathit{F}}_t^2 {\mathit{H}}_{f,2} \mathit{C}\left(\overline{n}\left({\alpha }_2\right)-{\alpha }_2\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_2\right)}{\mathit{D}}^T-{\overline{\lambda }}_1\overline{\mathit{X}}{\overline{\mathit{V}}}^1{\overline{\mathit{H}}}_1{\overline{\mathit{C}}}_1{\overline{\mathit{B}}}^1\overline{\mathit{D}}-{\overline{\lambda }}_2\overline{\mathit{X}}{\overline{\mathit{V}}}^2{\overline{\mathit{H}}}_2{\overline{C}}_2{\overline{\mathit{B}}}^2\overline{\mathit{D}}]\mathit{A}=\mathit{G}\hfill \end{array} \end{gathered}$$

We choose if $\overline{N}=5$, the approximate solution of the problem for $N=4, 10, 21,$ respectively

$$\begin{gathered} u_4\left(x\right)=-0.0569736258x^4-0.1665069441x^3-0.49426610634x^2-1.0017991281x^{ } \\ +0.001799128086 \end{gathered}$$

$$\begin{gathered} u_{10}\left(x\right)=-0.0112791176x^{10}+0.04361473251x^9-0.07167515469x^8+0.06498524957x^7 \\ \begin{array}{cc}\hfill & -0.03711020557x^6+0.00378578712x^5-0.0441096211x^4\hfill \\ \hfill & -0.1664357675x^3-0.4999839118x^2-0.999999027x^{ } \hfill \\ \hfill & -9.72893359856624\times {10}^{-7},\hfill \end{array} \end{gathered}$$

and

$$\begin{gathered} u_{21}\left(x\right)=2649.390397x^{21}-26487.723779x^{20}+123552.0169x^{19}-357222.928001x^{18} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}\hfill & +717359.135397x^{17}-1062540.52114x^{16}+1203203.42788x^{15} \hfill \\ \hfill & -1065423.12951x^{14}+748326.26838x^{13}-420462.50002x^{12} \hfill \\ \hfill & +189747.93897x^{11}-68792.817648x^{10}+19970.45119x^9-4609.812041x^8 \hfill \\ \hfill & +836.75503981x^7-117.510711061x^6+12.4665060412x^5 \hfill \\ \hfill & -1.0098682157x^4-0.11440971986x^3-0.50180747163x^2 \hfill \\ \hfill & -0.99996814632x^{ }-3.1859671582\times {10}^{-5}\hfill \end{array} \end{gathered}$$

we choose if $\overline{N}=10,$ the approximate solution of the problem for $N=4, 10, 21,$ respectively

$$u_4\left(x\right)=-0.05759757087x^4-0.1652553696x^3-0.4950425756x^2-1.001653973x^{ }+0.001653973$$

$$\begin{gathered} u_{10}\left(x\right)=3.399204141e-6x^{10}-2.0679446987e-5x^9+1.254303583e-5x^8 \\ \begin{array}{cc}\hfill & -0.000242474208x^7-0.001356395456x^6-0.008348879766x^5\hfill \\ \hfill & -0.0416618228x^4-0.166667629x^3-0.4999998846x^2\hfill \\ \hfill & -1.0000000056x+5.656145543\times {10}^{-9},\hfill \end{array} \end{gathered}$$

and

$$\begin{gathered} u_{21}\left(x\right)=0.0366737297x^{21}-0.36713669x^{20}+1.714854990x^{19}-4.9651991545x^{18} \\ \begin{array}{cc}\hfill & +9.9858420757x^{17}-14.8143137167x^{16}+16.8038086x^{15}\hfill \\ \hfill & -14.90654178x^{14}+10.49049965x^{13}-5.9069139812x^{12}\hfill \\ \hfill & +2.6719828271x^{11}-0.97127108855x^{10}+0.282792909x^9\hfill \\ \hfill & -0.0655249586x^8+0.0117379624x^7-0.00307302784x^6\hfill \\ \hfill & -0.008153528584x^5-0.041680721487x^4-0.1666659001x^3\hfill \\ \hfill & -0.5000000273x^2-0.9999999995x -4.670435338\times {10}^{-10}\hfill \end{array} \end{gathered}$$

Similarly doing it for $\overline{N}=16$, the approximate solution of the problem for $N=4, 10, 21,$ respectively

$$\begin{gathered} u_4\left(x\right)=-0.057597577x^4-0.165255357x^3-0.495042583x^2-1.001653971x^{ } \\ +0.0016539713 \end{gathered}$$

$$\begin{gathered} u_{10}\left(x\right)=3.4214653141e-6x^{10}-2.0612514452e-5x^9+1.2077128326e-5x^8 \\ \begin{array}{cc}\hfill & -0.000241586656x^7-0.0013572657x^6-0.0083483752x^5\hfill \\ \hfill & -0.0416620032x^4-0.1666675894x^3-0.499999889x^2\hfill \\ \hfill & -1.000000005x+5.395456491\times {10}^{-9}\hfill \end{array} \end{gathered}$$

and

$$\begin{gathered} u_{21}\left(x\right)=-0.000949182285x^{21}+0.00893091688x^{20}-0.03880287140x^{19} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}\hfill & +0.103101802x^{18}-0.186859983473284x^{17}+0.24354820561x^{16}\hfill \\ \hfill & -0.2337645213x^{15}+0.1652048762x^{14}-0.08287855594x^{13}\hfill \\ \hfill & +0.0253870484x^{12}-0.0005131248x^{11}-0.0043926365x^{10}\hfill \\ \hfill & +0.00281921684x^9-0.0010956089x^8+8.623002860e-5x^7\hfill \\ \hfill & -0.00144403537x^6-0.00832553568x^5-0.0416674549x^4-0.16666661x^3\hfill \\ \hfill & -0.500000002x^2-0.99999999995x^{ }-4.983565034\times {10}^{-11}.\hfill \end{array} \end{gathered}$$

In

Table 3. Comparison between the exact solution $u\left(x\right)$ and approximate solution $u_N\left(x\right)$ for Example 3.

(a)
$x_{\mathrm{𝒾}}$	Exact Solution Example 3.	$\overline{N}=5$
		$N$-Approximate Solution $u_N\left(x\right)$
		$N=4$	$N=10$	$N=21$
$0.0$	$0.0$	$0.00179913$	−9.72893360 × 10−7	$-3.1859672\times {10}^{-5}$
$0.1$	$-0.10517092$	$-0.10349565$	$-0.105171555$	$-0.10520278$
$0.2$	$-0.22140276$	$-0.21975456$	$-0.221402690$	$-0.22143482$
$0.3$	$-0.34985881$	$-0.34818173$	−0.349857839	$-0.34989091$
$0.4$	$-0.4918247$	$-0.49011807$	$-0.491822593$	$-0.49185648$
$0.5$	$-0.64872127$	$-0.64704118$	$-0.648717441$	$-0.64875193$
$0.6$	$-0.8221188$	$-0.82056543$	$-0.822111722$	$-0.82214657$
$0.7$	$-1.01375271$	$-1.0124419$	$-1.01373901$	−1.01377394
$0.8$	$-1.22554093$	$-1.22455843$	$-1.22551392$	$-1.22554873$
$0.9$	$-1.45960311$	$-1.45893959$	$-1.45955258$	$-1.45958517$
$1.0$	$-1.71828183$	$-1.71774668$	$-1.71820801$	$-1.71823455$
L.S.E.		$2.3130923\times {10}^{-05}$	$8.99107006\times {10}^{-09}$	$9.87875607\times {10}^{-09}$
Running Time/Sec		$0.171678066$	$0.49988222$	$2.92118477$
(b)
$x_{\mathrm{𝒾}}$	Exact Solution Example 3.	$\overline{N}=10$
		$N$-Approximate Solution $u_N\left(x\right)$
		$N=4$	$N=10$	$N=21$
$0$	$0.0$	$1.65397131\times {10}^{-03}$	$5.65614554\times {10}^{-09}$	$-4.67043534\times {10}^{-10}$
$0.1$	$-0.10517092$	$-1.03632867\times {10}^{-01}$	$-1.05170912\times {10}^{-01}$	$-1.05170919\times {10}^{-01}$
$0.2$	$-0.22140276$	$-2.19892725\times {10}^{-01}$	$-2.21402752\times {10}^{-01}$	−2.21402758 × 10−01
$0.3$	$-0.34985881$	$-3.48324488\times {10}^{-01}$	−3.49858802 × 10−01	$-3.49858808\times {10}^{-01}$
$0.4$	$-0.4918247$	$-4.90265271\times {10}^{-01}$	$-4.91824692\times {10}^{-01}$	$-4.91824698\times {10}^{-01}$
$0.5$	$-0.64872127$	$-6.47190428\times {10}^{-01}$	$-6.48721264\times {10}^{-01}$	$-6.48721271\times {10}^{-01}$
$0.6$	$-0.8221188$	−8.20713545 × 10−01	$-8.22118794\times {10}^{-01}$	$-8.22118801\times {10}^{-01}$
$0.7$	$-1.01375271$	$-1.01258644$	$-1.01375270$	$-1.01375271$
$0.8$	$-1.22554093$	$-1.22469917$	$-1.22554092$	$-1.22554093$
$0.9$	$-1.45960311$	$-1.45908002$	−1.45960311	$-1.45960311$
$1$	$-1.71828183$	$-1.71789552$	$-1.71828182$	$-1.71828183$
L.S.E.		$1.89772271\times {10}^{-05}$	$3.87412311\times {10}^{-16}$	$2.32392057\times {10}^{-18}$
Running Time/Sec		$0.187295436$	$0.515326499$	$2.749377012$
(c)
$x_{\mathrm{𝒾}}$	Exact Solution Example 3.	$\overline{N}=16$
		$N$-Approximate Solution $u_N\left(x\right)$
		$N=4$	$N=10$	$N=21$
$0.0$	$0.$	1.65397131 × 10−03	$5.39545649\times {10}^{-09}$	$-4.98356503\times {10}^{-11}$
$0.1$	$-0.10517092$	$-1.03632867\times {10}^{-01}$	$-1.05170913\times {10}^{-01}$	$-1.05170918\times {10}^{-01}$
$0.2$	$-0.22140276$	$-2.19892725\times {10}^{-01}$	$-2.21402753\times {10}^{-01}$	$-2.21402758\times {10}^{-01}$
$0.3$	$-0.34985881$	$-3.48324488\times {10}^{-01}$	$-3.49858802\times {10}^{-01}$	$-3.49858808\times {10}^{-01}$
$0.4$	$-0.4918247$	$-4.90265271\times {10}^{-01}$	$-4.91824692\times {10}^{-01}$	$-4.91824698\times {10}^{-01}$
$0.5$	$-0.64872127$	$-6.47190428\times {10}^{-01}$	$-6.48721265\times {10}^{-01}$	$-6.48721271\times {10}^{-01}$
$0.6$	$-0.8221188$	$-8.20713545\times {10}^{-01}$	$-8.22118794\times {10}^{-01}$	$-8.22118800\times {10}^{-01}$
$0.7$	$-1.01375271$	$-1.01258644$	$-1.01375270$	$-1.01375271$
$0.8$	$-1.22554093$	$-1.22469917$	−1.22554092	$-1.22554093$
$0.9$	$-1.45960311$	$-1.45908002$	$-1.45960311$	$-1.45960311$
$1.0$	$-1.71828183$	$-1.71789552$	$-1.71828182$	$-1.71828183$
L.S.E.		$1.89771908\times {10}^{-05}$	$3.4403856\times {10}^{-16}$	$3.45079782\times {10}^{-20}$
Running Time/Sec		$0.203105688$	$0.5199816226$	$2.90862059$

presents a comparison between the exact solution $u\left(x\right)$ and approximate solution $u_N\left(x\right)$, when we choose $\overline{N}=5,10,\mathrm{and}16,$ respectively. For each of them we chose $N=4,10,\mathrm{and}21,$ respectively depending on the least square error and running time.

Figure 3a–c illustrates a comparison between the exact solution and approximate solution of (LF-VIFDEs) of equation above, respectively. To show the result of the proposed method to an exact solution, we present Table 3, respectively. Each of the plots is drawn with our Python program version 3.8.8 (2021).

Example 4.

Suppose that the following linear Fredholm–Volterra fractional integro-differential equation given by

$$\begin{gathered} {}_0{}^{C}D{}_x^{0.8}u\left(x\right)=g\left(x\right)+\underset{0}{\overset{2}{{\displaystyle \int }}}\left(\frac{1}{2}\left(t-2\right){}_0{}^{C}D{}_t^{0.3}u\left(t\right)+\left(t+\cos \left(x\right)\right){}_0{}^{C}D{}_t^{1.3}u\left(t\right)\right)dt \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\underset{0}{\overset{x}{{\displaystyle \int }}}\left(\frac{1}{4}\left(\frac{x-t}{2}\right)u\left(t\right)+\left[\tan \left(x\right)t\right]{}_0{}^{C}D{}_t^{2.7}u\left(t\right)\right)dt, 0\le x,t\le 2 \end{gathered}$$

where

$$\begin{gathered} g\left(x\right)=\frac{2}{\mathsf{\Gamma }\left(2.2\right)}{x}^{1.2}-\frac{1}{\mathsf{\Gamma }\left(1.2\right)}{x}^{0.2}-\frac{\left(2.7\right)}{\mathsf{\Gamma }\left(4.7\right)}{2}^{3.7}-\frac{1}{\mathsf{\Gamma }\left(3.7\right)}\left(\left(1.7\right)2^{3.7}-2^{3.7}-\left(1.7\right)2^{1.7}\right) \\ -\frac{1}{\mathsf{\Gamma }\left(2.7\right)}\left(2^{1.7}+2^{2.7}\cos \left(x\right)\right)-\frac{1}{16}\left(\frac{1}{6}{x}^4-\frac{1}{3}{x}^3+x^2\right) \end{gathered}$$

with the boundary conditions

$$u^{ }\left(0\right)+u^{ }\left(2\right)=4, u^{\left(1\right)}\left(0\right)+u^{\left(1\right)}\left(2\right)=2 and u^{\left(2\right)}\left(0\right)+u^{\left(2\right)}\left(2\right)=4$$

which is the exact solution $u\left(x\right)=x^2-x+1.$

Now let us find the approximate solution given by the $N$-truncated Bessel series

$$u\left(x\right)\cong u_N\left(x\right)={\displaystyle \sum }_{𝓇=0}^N{a}_{𝓇}{J}_{𝓇}\left(x\right)$$

Here, from consider example we have:

$$\begin{array}{cc}{\sigma }_1=0.8\to \overline{n}\left({\sigma }_1\right)={\sigma }_1=1\hfill & \hfill \\ {\alpha }_1=0.3,{\alpha }_2=1.3\to \overline{n}\left({\alpha }_1\right)={\alpha }_1=1,\overline{n}\left({\alpha }_2\right)={\alpha }_2=2\hfill & \hfill \\ {\beta }_0=0,{\beta }_1=2.7 \to \overline{n}\left({\beta }_0\right)={\beta }_0=0, \overline{n}\left({\beta }_1\right)={\beta }_1=3\hfill & \hfill \\ \mu =\max \left\{\lceil 0.8\rceil ,\lceil 1.3\rceil ,\lceil 2.7\rceil \right\}=3 ,{\lambda }_1=\frac{1}{2},{\lambda }_2=1, {\overline{\lambda }}_0=\frac{1}{4}, {\overline{\lambda }}_1=1\hfill \end{array}$$

and ${\mathit{p}}_1\left(x\right)=0,F_1\left(x,t\right)=t-2,F_2\left(x,t\right)=t+\cos \left(x\right),V_0\left(x,t\right)=\left(\frac{x-t}{2}\right),V_1\left(x,t\right)=\tan \left(x\right)t,$ from Equation (31), the fundamental matrix equation of the given problem is written as

$$\begin{gathered} [\mathit{y}\left(1\right) \mathit{X} \mathit{C}\left(\overline{n}\left({\sigma }_1\right)-{\sigma }_1\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\sigma }_1\right)}{\mathit{D}}^T-{\lambda }_1\mathit{X}{\mathit{F}}_t^1 {\mathit{H}}_{f,1} \mathit{C}\left(\overline{n}\left({\alpha }_1\right)-{\alpha }_1\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_1\right)}{\mathit{D}}^T \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\lambda }_2\mathit{X}{\mathit{F}}_t^2 {\mathit{H}}_{f,2} \mathit{C}\left(\overline{n}\left({\alpha }_2\right)-{\alpha }_2\right) {\left({\mathit{B}}^T\right)}^{\overline{n}\left({\alpha }_2\right)}{\mathit{D}}^T-{\overline{\lambda }}_0\overline{\mathit{X}}{\overline{\mathit{V}}}^0{\overline{\mathit{H}}}_0\overline{\mathit{D}}-{\overline{\lambda }}_1\overline{\mathit{X}}{\overline{\mathit{V}}}^1{\overline{\mathit{H}}}_1{\overline{\mathit{C}}}_1{\overline{\mathit{B}}}^1\overline{\mathit{D}}]\mathit{A}=\mathit{G} \end{gathered}$$

Thus, the approximate solution of the problem for $N=5, 9, 12,$ respectively

$$\begin{gathered} u_5\left(x\right)=0.00595039974990515x^5-0.0295085156834566x^4+0.0436273982107945x^3 \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+0.98520400357289x^2-0.998052135471447x^{ }+0.99399626495166, \end{gathered}$$

$$\begin{gathered} u_9\left(x\right)=-6.93792068979952\times {10}^{-7}{x}^9+6.12616976553793\times {10}^{-6}{x}^8-2.13390478665744\times {10}^{-5}{x}^7 \\ \begin{array}{cc}\hfill & +3.69079850914571\times {10}^{-5}{x}^6-3.01713615634192\times {10}^{-5}{x}^5\hfill \\ \hfill & +1.16291867668927\times {10}^{-6}{x}^4+1.66943823538962\times {10}^{-5}{x}^3\hfill \\ \hfill & +0.999993226980557x^2-0.999998942794685x^{ }+0.99999725431996,\hfill \end{array} \end{gathered}$$

and

$$\begin{gathered} u_{12}\left(x\right)=-5.96684668513465\times {10}^{-10}{x}^{12}+6.51569719257085\times {10}^{-9}{x}^{11} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{cc}\hfill & -3.10471674650442\times {10}^{-8}{x}^{10}+8.49087433427835\times {10}^{-8}{x}^9\hfill \\ \hfill & -1.47156926969634\times {10}^{-7}{x}^8+1.67915537193858\times {10}^{-7}{x}^7\hfill \\ \hfill & -1.26723286902409\times {10}^{-7}{x}^6+6.23462096628266\times {10}^{-8}{x}^5\hfill \\ \hfill & -2.05277535525461\times {10}^{-8}{x}^4+5.82167120066757\times {10}^{-9}{x}^3\hfill \\ \hfill & +0.999999998891922x^2-0.999999999885384x^{ }+0.999999999582705.\hfill \end{array} \end{gathered}$$

In

Table 4. Comparison between the exact solution $u\left(x\right)$ and approximate solution $u_N\left(x\right)$ for Example 4.

${\mathit{x}}_{\mathbf{𝒾}}$	Exact Solution Example 4.	$\mathit{N}$$\text{-}\mathbf{Approximate}\mathbf{Solution}{\mathit{u}}_{\mathit{N}}\left(\mathit{x}\right)$
		$\mathit{N}=\mathbf{5}$	$\mathit{N}=\mathbf{9}$	$\mathit{N}=\mathbf{12}$
$0.0$	$1.00$	$0.99399626$	$0.99999725$	$1.00$
$0.1$	$0.91$	$0.90408383$	$0.90999731$	$0.91$
$0.2$	$0.84$	$0.83409771$	$0.83999732$	$0.84$
$0.3$	$0.79$	$0.78420236$	$0.78999737$	$0.79$
$0.4$	$0.76$	$0.75450572$	$0.7599975$	$0.76$
$0.5$	$0.75$	$0.74506629$	$0.74999774$	$0.75$
$0.6$	$0.76$	$0.75590034$	$0.75999808$	$0.76$
$0.7$	$0.79$	$0.78698902$	$0.78999852$	$0.79$
$0.8$	$0.84$	$0.83828549$	$0.83999904$	$0.84$
$0.9$	$0.91$	$0.90972207$	$0.90999961$	$0.91$
$1.0$	$1.00$	$1.00121742$	$1.00000023$	$1.00$
L.S.E.		$0.0002244$	$4.72053254\times {10}^{-11}$	$1.11096853\times {10}^{-18}$
Running Time/Sec		$0.568155527$	$1.649296999$	$6.3344522$

. presents a comparison between the exact solution $u\left(x\right)$ and approximate solution $u_N\left(x\right)$ for $N=5,9and12,$ respectively, depending on the least square error and running time.

Figure 4 illustrates a comparison between the exact solution and approximate solution of linear (FVIFDEs). To show the result of the proposed method to an exact solution, we present Table 4. Each of the plots is drawn with our Python program version 3.8.8 (2021).

6. Conclusions

Multi-fractional order linear integro-differential equations are generally difficult to solve analytically. In many situations, it is necessary to approximate solutions. In this work, we present a new technique for numerically solving the linear Fredholm–Volterra integro-fractional differential equation of multi-fractional order of the Caputo sense using first-order Bessel polynomials. The comparison of the results achieved with the exact solution, the exact solution, and the other methods suggests that the procedure is very effective and convenient. We introduced this with some illustrative examples of the approach and their least square error to minimize the error terms on the specified domain and running time are also given in tabular form. It is obvious that as $N$ rises, the error rate reduces and the answer becomes closer to the exact solution. One significant benefit of the technique is that the Bessel coefficients of the solution may be determined relatively quickly using computer code developed in Python v3.8.8 (2021). As an example, consider the Python v3.8.8 (2021).

Future directions: Using the residual error function, we can enhance the Bessel collocation method for solving the multi-high fractional-order system of Fredholm–Volterra integro-differential equations and their delay. This technique can also be used to make an accurate error estimation.