**1. Introduction**

Fractional derivatives and integrals have recently gained high interest in many fields of science. The ability of classifying and capturing the memory and hereditary properties of various materials and processes is an advantage of fractional derivatives over their integer counterparts, e.g., the modeling of anomalous diffusion by fractional differential equations gives more informative and interesting models [1]. For time-fractional differential equations, the memory feature implies that all previous information is needed to evaluate the time fractional derivative at the current time level. Accordingly, designing a numerical differentiation formula of good accuracy is as ever paramount, but especially hard. The approximation formulas based on the interpolation approximation, such as *L*1 [2] and *L*<sup>2</sup> − 1*<sup>σ</sup>* [3], are of significance to design numerical algorithms to solve time-fractional differential equations. Demonstrated applications in numerous seemingly diverse and widespread fields of physics, such as in porous and glassy materials, in percolation clusters over fractals to semi-conductors, polymers, random media, and beyond, like geophysical and biological systems or processes (e.g., [4,5]), can be effectively modeled by time-fractional diffusion-wave equations of different types. Here, we are seeking to design a compact difference scheme that behaves linearly to numerically solve the non-linear

delayed multiterm time fractional convection diffusion-wave equation (dmfCDWEs) with spatial variable coefficients. More specifically, we consider

$$\sum\_{r=0}^{m} p\_r \frac{\partial^{\mu\_t} u(\mathbf{x}, t)}{\partial t^{\mu\_t}} = \frac{\partial}{\partial \mathbf{x}} \left( q\_1(\mathbf{x}) \frac{\partial u}{\partial \mathbf{x}} \right) + q\_2(\mathbf{x}) \frac{\partial u}{\partial \mathbf{x}} + f(u(\mathbf{x}, t), u(\mathbf{x}, t-\mathbf{s}), \mathbf{x}, t), \ 0 < t \le T, \ 0 \le \mathbf{x} \le L,\tag{1a}$$

with the following initial and boundary conditions

$$u(\mathbf{x},t) = d(\mathbf{x},t), \quad 0 \le \mathbf{x} \le L, \quad t \in [-s,0), \quad \frac{\partial u(\mathbf{x},0)}{\partial t} = \boldsymbol{\upvarphi}(\mathbf{x}) = \lim\_{t \to -0} \frac{\partial d(\mathbf{x},t)}{\partial t}, \tag{1b}$$

$$u(0,t) = \phi\_0(t), \quad u(L,t) = \phi\_L(t), \quad 0 < t \le T,\tag{1c}$$

where *s* > 0 is a fixed delay parameter, *q*1(*x*) and *q*2(*x*) are functions chosen, respectively, to be sufficiently and arbitrary differentiable functions. The fractional derivative is defined in Caputo sense and the fractional orders {*α<sup>r</sup>* | 1 ≤ *r* ≤ *m*} are specified in the manner {1 < *α<sup>m</sup>* ≤ *αm*−<sup>1</sup> < ··· < *α*<sup>0</sup> = 2}. The existence and uniqueness of the global mild solutions for the problem of nonlinear fractional reaction–diffusion equations with delay and Caputo's fractional derivatives are addressed in [6].

This work can be considered to be an extension of our previously published work [7], in which we discussed a single term time fractional wave equation with spatial constant coefficients. The scheme was of 2 − *α* order in time and fourth in space. Here, we treat the multiterm time fractional order case with spatial variable coefficients and seek to have temporal second order of convergence. Accordingly, we hinged on the proposed numerical formula in [8] to approximate the multiterm Caputo fractional derivatives of order *αr*(0 < *α<sup>r</sup>* ≤ 1) at the super-convergent point. The formula *L*<sup>2</sup> − 1*<sup>σ</sup>* can achieve at least second-order accuracy at this point. We rely on the compact operator proposed in [9] in order to attain a fourth order accuracy with respect to space in case of having spatial variable coefficients.

There are some results and findings available regarding the theoretical analysis and numerical computation of single term time fractional sub or super diffusion equations with delay. In [10], the authors introduced a satisfactory numerical method for time fractional diffusion equations with delay. In [11], a novel discrete Grönwall inequality is used to simplify the analysis of difference schemes for time-fractional multi-delayed sub-diffusion equations. Convergence and stability of a compact finite difference method for nonlinear time fractional reaction-diffusion equations with a fixed delay are proposed in [12] by the aid of a new discrete form of the fractional Grönwall inequalities. A numerical solution for a class of time fractional diffusion equations with delay is proposed in [13] that is based on a smooth difference approximation of specific *L*1 type. Additionally, there are many difference and spectral approaches proposed for multiterm or distributed order time fractional differential equations. An efficient spectral method that is based on Jacobi–Gauss–Radau collocation is applied in order to solve a system of multi-dimensional distributed-order generalized Schrödinger equations in [14]. A combined difference and Galerkin–Legendre spectral method in [15] is used to solve time fractional diffusion equations with nonlinear source term. A Legendre spectral-collocation method for the numerical solution of distributed-order fractional initial value problems is designed in [16]. In [17], the authors proposed a spectral *τ*-scheme to discretize the fractional diffusion equation with distributed-order fractional derivative in time and Dirichlet boundary conditions. The model solution is expanded in multi-dimensions in terms of Legendre polynomials and the discrete equations are obtained with the *τ*-method. The two-dimensional distributed-order time fractional cable equation is numerically solved based on the finite difference/spectral method, as clarified in [18]. Two classes of finite difference methods that are based on backward differential formula discretization in the temporal direction are proposed in [19] to efficiently solve the semilinear space fractional reaction–diffusion equation with time-delay. The coefficients of the problem are constants, a fractional centered difference approximation is employed for the space fractional derivative [20], and it gains a fourth order approximation in space due to the use of a specific compact operator [21].

The main purpose of our work is the manufacturing of a difference scheme for problems of the kind of (1). Until now, few works hace paid close attention to the multiterm fractional wave equation with variable coefficients and delay argument simultaneously. It is well known that it may be a more challenging task to solve the fractional partial differential equation with delay effectively, since the evolution of a fractional partial differential equation with delay at time *t* not only depends on its value at *t* − *s*, but also depends on all previous solutions due to the non-locality of the fractional operator. Higher order numerical schemes are extremely scarce and difficult in regard to the analysis and implementation for the variable coefficient multiterm time fractional convection-diffusion wave equation with delay. For a single temporal term for that kind of problem, we refer to [22]. For the multiterm problem, we design the difference scheme at a super-convergent point to gain a high order of convergence [8] which is more challenging than the single term problem at the level of theoretical analysis.

In order to simplify the problem, an exponential transformation technique [23] can be applied to system (1) to avoid the difficulties resulting from the spatial variable when constructing high order compact difference methods. That technique is used to eliminate the convection term. Assume that *q*1(*x*) *<sup>q</sup>*2(*x*) is integrable in the spatial interval [0, *<sup>L</sup>*] and let *<sup>u</sup>*(*x*, *<sup>t</sup>*) = exp *x* <sup>0</sup> *r*˜(**s**)*d***s** *ω*(*x*, *t*), this yields

$$\frac{\partial^{a\_r}u}{\partial t^{a\_r}} = \exp\left(\int\_0^x \overline{r}(\mathbf{s})d\mathbf{s}\right)\frac{\partial^{a\_r}\omega}{\partial t^{a\_r}},\tag{2}$$

$$\frac{\partial \mu}{\partial \mathbf{x}} = \exp\left(\int\_0^\mathbf{x} \vec{r}(\mathbf{s}) d\mathbf{s}\right) \left[\vec{r}(\mathbf{x})\omega(\mathbf{x}, t) + \frac{\partial \omega}{\partial \mathbf{x}}\right],\tag{3}$$

$$\frac{\partial^2 \underline{u}}{\partial \mathbf{x}^2} = \exp\left(\int\_0^\mathbf{x} \overline{r}(\mathbf{s}) d\mathbf{s}\right) \left[\overline{r}^2(\mathbf{x})\omega(\mathbf{x},t) + 2\overline{r}(\mathbf{x})\frac{\partial \omega}{\partial \mathbf{x}} + \overline{r}'(\mathbf{x})\omega(\mathbf{x},t) + \frac{\partial^2 \omega}{\partial \mathbf{x}^2}\right].\tag{4}$$

When substituting (2)–(4) into (1), one directly obtains

$$\begin{split} \sum\_{r=0}^{m} p\_r \frac{\partial^{\mathbf{z}\_r} \omega(\mathbf{x}, t)}{\partial t^{\mathbf{z}\_r}} &= q\_1(\mathbf{x}) \frac{\partial^2 \omega}{\partial \mathbf{x}^2} + \left[ q\_1'(\mathbf{x}) + 2q\_1(\mathbf{x}) \widetilde{\mathbf{r}}(\mathbf{x}) + q\_2(\mathbf{x}) \right] \frac{\partial \omega}{\partial \mathbf{x}} \\ &+ \left[ q\_1(\mathbf{x}) \widetilde{\mathbf{r}}^2(\mathbf{x}) + q\_1(\mathbf{x}) \widetilde{\mathbf{r}}'(\mathbf{x}) + q\_1'(\mathbf{x}) \widetilde{\mathbf{r}}(\mathbf{x}) + q\_2(\mathbf{x}) \widetilde{\mathbf{r}}(\mathbf{x}) \right] \omega(\mathbf{x}, t) \\ &+ \exp\left( -\int\_0^\mathbf{x} \widetilde{\tau}(\mathbf{s}) d\mathbf{s} \right) \mathbf{g}(\omega(\mathbf{x}, t), \omega(\mathbf{x}, t-\mathbf{s}), \mathbf{x}, t), \end{split}$$

where, *<sup>g</sup>*(*ω*(*x*, *<sup>t</sup>*), *<sup>ω</sup>*(*x*, *<sup>t</sup>* <sup>−</sup> *<sup>s</sup>*), *<sup>x</sup>*, *<sup>t</sup>*) = *<sup>f</sup>*(*u*(*x*, *<sup>t</sup>*), *<sup>u</sup>*(*x*, *<sup>t</sup>* <sup>−</sup> *<sup>s</sup>*), *<sup>x</sup>*, *<sup>t</sup>*). By selecting *<sup>r</sup>*˜(*x*) = <sup>−</sup><sup>1</sup> 2 *q*2(*x*) *<sup>q</sup>*1(*x*), the system (1) is transformed in

$$\sum\_{r=0}^{m} p\_r \frac{\partial^{\omega\_r} \omega(\mathbf{x}, t)}{\partial t^{\omega\_r}} = \frac{\partial}{\partial \mathbf{x}} \left( q\_1(\mathbf{x}) \frac{\partial \omega}{\partial \mathbf{x}} \right) + \tilde{f}(\omega(\mathbf{x}, t), \omega(\mathbf{x}, t-\mathbf{s}), \mathbf{x}, t), \ 0 < t \le T, \ 0 \le \mathbf{x} \le L,\tag{5a}$$

with the following initial and boundary conditions

$$\omega(\mathbf{x},t) = \exp\left(-\int\_0^\mathbf{x} \mathfrak{f}(\mathbf{s})d\mathbf{s}\right)d(\mathbf{x},t) := \tilde{d}(\mathbf{x},t), \quad t \in [-s,0), \quad \frac{\partial \omega(\mathbf{x},0)}{\partial t} = \lim\_{t \to -0} \frac{\partial \tilde{d}(\mathbf{x},t)}{\partial t} := \tilde{\mathfrak{y}}(\mathbf{x}), \tag{5b}$$

$$
\omega(0,t) = \phi\_0(t), \quad \omega(L,t) = \exp\left(-\int\_0^T \mathbb{P}(\mathbf{s})d\mathbf{s}\right)\phi\_L(t) := \phi\_L(t), \quad 0 < t \le T,\tag{5c}
$$

where

$$\begin{split} f(\omega(\mathbf{x},t),\omega(\mathbf{x},t-\mathbf{s}),\mathbf{x},t) &= \left[q\_1(\mathbf{x})\overline{r}^2 + q\_1(\mathbf{x})\overline{r}(\mathbf{x}) + q\_1'(\mathbf{x})\overline{r}(\mathbf{x}) + q\_2(\mathbf{x})\overline{r}(\mathbf{x})\right] \omega(\mathbf{x},t) \\ &+ \exp\left(-\int\_0^\mathbf{x} \overline{r}(\mathbf{s})d\mathbf{s}\right) \mathbf{g}(\omega(\mathbf{x},t),\omega(\mathbf{x},t-\mathbf{s}),\mathbf{x},t). \end{split}$$

In order to transform (5) to a system with zero Dirichlet boundary conditions, we define *h*(*x*, *t*) := *φ*0(*t*) + *<sup>x</sup> <sup>L</sup>* (*φ*˜*L*(*t*) − *<sup>φ</sup>*0(*t*)) and introduce the new function *<sup>ν</sup>*(*x*, *<sup>t</sup>*) = *<sup>ω</sup>*(*x*, *<sup>t</sup>*) − *<sup>h</sup>*(*x*, *<sup>t</sup>*). Hence, we have

$$\sum\_{r=0}^{m} p\_r \frac{\partial^{a\_r} \nu(\mathbf{x}, t)}{\partial t^{a\_r}} = \frac{\partial}{\partial \mathbf{x}} \left( q\_1(\mathbf{x}) \frac{\partial \nu}{\partial \mathbf{x}} \right) + \tilde{f}(\nu(\mathbf{x}, t), \nu(\mathbf{x}, t-\mathbf{s}), \mathbf{x}, t), \ 0 < t \le T, \ 0 \le \mathbf{x} \le L,\tag{6a}$$

with the following initial and boundary conditions

$$w(\mathbf{x},t) = \mathbf{f}(\mathbf{x},t), \quad 0 \le \mathbf{x} \le L, \quad t \in [-s,0), \quad \frac{\partial v(\mathbf{x},0)}{\partial t} = \hat{\psi}(\mathbf{x}) = \lim\_{t \to -0} \frac{\partial r(\mathbf{x},t)}{\partial t}, \tag{6b}$$

$$
\upsilon(0,t) = \upsilon(L,t) = 0, \quad t > 0. \tag{6c}
$$

According to the new transformed system (6), we analytically overcame the first degree of complexity by the elimination of the convection term. There are still two degrees of complexity to be numerically overcome; the multiterm fractional order on the one hand and the nonlinear delay on the other.

Throughout this work, we assume that the function *f*(*μ*, *v*, *x*, *t*) and the solution *ν*(*x*, *t*) of (6) are sufficiently smooth in the following sense:


$$\mathcal{L}\_1 = \sup\_{0 < \mathbf{x} < L, \mathbf{0} < t \le T} \left| \tilde{f}\_{\mu}(\nu(\mathbf{x}, t) + \boldsymbol{\varepsilon}\_1, \nu(\mathbf{x}, t - s) + \boldsymbol{\varepsilon}\_2, \mathbf{x}, t) \right|, \tag{7a}$$

$$\begin{array}{c} |c\_1| \le \varepsilon\_0, |c\_2| \le \varepsilon\_0\\ \varepsilon\_2 = \max\_{\substack{0 < \mathbf{x} < L, \, 0 < t \le T\\ |\mathbf{c}\_1| \le \varepsilon\_0, |\mathbf{c}\_2| \le \varepsilon\_0}} \left| f\_\nu(u(\mathbf{x}, t) + \varepsilon\_1, \nu(\mathbf{x}, t - \mathbf{s}) + \varepsilon\_2, \mathbf{x}, t) \right|. \end{array} \tag{7b}$$

The structure of this paper is arranged, as follows. First, we introduce a derivation of the compact difference scheme. Next, in the third section, convergence and stability for the compact difference scheme are carried out. Finally, the paper ends with a numerical illustration and a conclusion.
