*2.1. L*<sup>2</sup> − 1*<sup>σ</sup> Super-Convergence Scheme*

Here, we give a preliminary for the Alikhanov scheme. Denote *γ<sup>r</sup>* = *α<sup>r</sup>* − 1 (0 ≤ *r* ≤ *m*) and

$$\mathcal{F}(\sigma) = \sum\_{r=0}^{m} \frac{p\_r}{\Gamma(3-\gamma\_r)} \sigma^{1-\gamma\_r} \left[\sigma - (1-\frac{\gamma\_r}{2})\right] \pi^{2-\gamma\_r}.$$

such that 0 < *γ<sup>m</sup>* < *γm*−<sup>1</sup> < ··· < *γ*<sup>0</sup> ≤ 1. Additionally, now we invoke the following lemmas from Alikhanov work.

**Lemma 2.1** ([8])**.** *A unique positive root σ*<sup>∗</sup> ∈ [1 − *γ*0/2, 1 − *γm*/2] *exists for* F(*σ*) = 0*.*

**Lemma 2.2** ([8])**.** *For m* = 0*, the root of* F(*σ*) = 0 *is* 1 − *γ*0/2*. However, if m* ≥ 1*, the root σ*<sup>∗</sup> *can be obtained by the Newton iteration method. The Newton iteration sequence* {*σk*}<sup>∞</sup> *<sup>k</sup>*=<sup>0</sup> *generated by σ*<sup>0</sup> = 1 − *γm*/2 *and <sup>σ</sup>k*+<sup>1</sup> <sup>=</sup> *<sup>σ</sup><sup>k</sup>* <sup>−</sup> <sup>F</sup>(*σ<sup>k</sup>* ) <sup>F</sup>(*σ<sup>k</sup>* ) *for k* <sup>=</sup> 0, 1, 2, . . . *is monotonically decreasing and convergent to <sup>σ</sup>*∗*.*

For the sake of simplicity, let *σ* = *σ*∗ here and later. Define [3]

$$a\_0 = \sigma^{1-\gamma}, \quad a\_l^\gamma = (\sigma+l)^{1-\gamma} - (\sigma+l-1)^{1-\gamma},$$

$$b\_l^\gamma = \frac{1}{2-\gamma} \left[ (l+\sigma)^{2-\gamma} - (l-1+\sigma)^{2-\gamma} \right] - \frac{1}{2} \left[ (l+\sigma)^{1-\gamma} + (l-1+\sigma)^{1-\gamma} \right],$$

for each *<sup>l</sup>* <sup>∈</sup> <sup>N</sup>+. Next, we define {*C*(*k*+1,*γ*) *<sup>n</sup>* }, as follows

$$\mathbb{C}\_{n}^{(k+1,\gamma)} = \begin{cases} a\_{0\prime} & n=0, \ k=0, \\ a\_{0}^{\gamma} + b\_{1\prime}^{\gamma} & n=0, \ k \ge 1, \\ a\_{n}^{\gamma} + b\_{n+1}^{\gamma} - b\_{n}^{\gamma} & 1 \le n \le k-1, \ k \ge 1, \\ a\_{k}^{\gamma} - b\_{k}^{\gamma} & n=k, \ k \ge 1. \end{cases} \tag{8}$$

Denote

$$\hat{\mathsf{C}}\_{n}^{(k+1)} = \sum\_{r=0}^{m} p\_r \frac{\mathsf{r}^{-\gamma\_r}}{\Gamma(2-\gamma\_r)} \mathsf{C}\_n^{(k+1,\gamma\_r)}, \quad \hat{b}\_n = \sum\_{r=0}^{m} p\_r \frac{\mathsf{r}^{-\gamma\_r}}{\Gamma(2-\gamma\_r)} b\_n^{(\gamma\_r)}, \quad n = 0, 1, \dots, k.$$

The next two lemmas are devoted to the properties of the coefficients *C*ˆ(*k*) *<sup>n</sup>* and ˆ *bn*.

**Lemma 2.3** ([8])**.** *Given any non-negative integer m and positive constants p*0, *p*1, ... , *pm, for any* {*γ<sup>r</sup>* ∈ (0, 1] | *r* = 0, 1, ··· , *m*} *it holds*

$$\mathfrak{C}\_1^{(k+1)} > \mathfrak{C}\_2^{(k+1)} > \dots > \mathfrak{C}\_{k-2}^{(k+1)} > \mathfrak{C}\_{k-1}^{(k+1)} > \sum\_{r=0}^m p\_r \frac{\mathfrak{r}^{-\gamma\_r}}{\Gamma(2-\gamma\_r)} \frac{1-\gamma\_r}{2} (k-1+\sigma)^{-\gamma\_r}.$$

*In addition, there exists a <sup>τ</sup>*<sup>0</sup> <sup>&</sup>gt; <sup>0</sup>*, such that* (2*<sup>σ</sup>* <sup>−</sup> <sup>1</sup>)*C*ˆ(*k*+1) <sup>0</sup> <sup>−</sup> *<sup>σ</sup>C*ˆ(*k*+1) <sup>1</sup> > 0, *when τ* ≤ *τ*0, *n* = 2, 3, ··· *, and C*ˆ(*k*+1) <sup>0</sup> <sup>&</sup>gt; *<sup>C</sup>*ˆ(*k*+1) <sup>1</sup> .

**Lemma 2.4** ([24])**.** *The sequences <sup>C</sup>*ˆ(*k*) *<sup>n</sup> and* ˆ *bn satisfy*

$$\mathcal{C}\_{n}^{(k+1)} = \begin{cases} \mathcal{C}\_{n}^{(k)} \, & 0 \le n \le k - 2, \\\mathcal{C}\_{n}^{(k)} + \hat{b}\_{n+1 \prime} \, & n = k - 1. \end{cases} \tag{9}$$

*Additionally, the following estimates hold:*

$$\sum\_{n=1}^{k} \mathbb{C}\_{n}^{(k+1)} \le \sum\_{r=0}^{m} p\_r \frac{3\pi^{-\gamma\_r}}{2\Gamma(2-\gamma\_r)} (k+\sigma)^{1-\gamma\_r}$$

*and*

$$\sum\_{m=1}^k \hat{b}\_m \le \sum\_{r=0}^m p\_r \frac{\gamma\_r \pi^{-\gamma\_r}}{2\Gamma(3-\gamma\_r)} (k+\sigma)^{1-\gamma\_r}.$$

Let W*<sup>h</sup>* = \$ *<sup>w</sup>* : <sup>Ω</sup>*h<sup>τ</sup>* <sup>→</sup> <sup>R</sup> <sup>|</sup> *<sup>w</sup>*(*xi*, *tk*) = *<sup>w</sup><sup>k</sup> <sup>i</sup>* ; *i* = 0, 1, . . . , *M*; *k* = −*n*0, −*n*<sup>0</sup> + 1, . . . , *N* % be a grid function space on <sup>Ω</sup>*h<sup>τ</sup>* and also define <sup>W</sup><sup>ˆ</sup> *<sup>h</sup>* <sup>=</sup> {*<sup>w</sup>* ∈ W*h*, *<sup>w</sup>*(*x*0, ·) = *<sup>w</sup>*(*xM*, ·) = <sup>0</sup>}. Introduce the following notations

*Mathematics* **2020**, *8*, 1696

$$w\_{i+1/2}^k = \frac{1}{2} \left[ w\_{i+1}^k + w\_i^k \right], \quad \delta\_\mathbf{x} w\_{i+1/2}^k = \frac{1}{h} \left[ w\_{i+1}^k - w\_i^k \right], \quad \delta\_\mathbf{l} w\_i^\frac{1}{2} = \frac{1}{\tau} \left[ w\_i^1 - w\_i^0 \right], \tag{10}$$

$$
\delta\_{\hat{x}} w\_i^k = \frac{1}{2h} \left[ w\_{i+1}^k - w\_{i-1}^k \right], \\
\delta\_{\hat{x}}^2 w\_i^k = \frac{1}{h} \left[ \delta\_{\hat{x}} w\_{i+1/2}^k - \delta\_{\hat{x}} w\_{i-1/2}^k \right], \tag{11}
$$

$$w\_i^{k+\sigma} = \sigma w\_i^{k+1} + (1 - \sigma)w\_i^k, \quad \partial\_l w\_i^k = \frac{1}{2\pi} \left[ (2\sigma + 1)w\_i^{k+1} - 4\sigma w\_i^k + (2\sigma - 1)w\_i^{k-1} \right]. \tag{12}$$

Moreover, denote κ<sup>1</sup> = (*q* <sup>1</sup>)2/*q*<sup>1</sup> <sup>−</sup> <sup>1</sup> <sup>2</sup> *q* <sup>1</sup> and <sup>κ</sup><sup>2</sup> <sup>=</sup> *<sup>q</sup>*<sup>1</sup> <sup>−</sup> *<sup>h</sup>*<sup>2</sup> <sup>12</sup>κ1. The compact operator acting on the spatial variable is defined as [9],

$$\mathcal{A}w\_{i} = \begin{cases} w\_{i} + \frac{h^{2}}{12} \left[ \delta\_{\hat{x}}^{2} w\_{i} - \delta\_{\hat{x}} \left( \frac{q\_{1}^{\prime}}{q\_{1}} w \right)\_{i} \right], & i = 1, \dots, M - 1, \\ w\_{i\prime} & i = 0, M. \end{cases}$$

**Lemma 2.5** ([9])**.** *Let g*(*x*) <sup>∈</sup> *<sup>C</sup>*6[0, *<sup>L</sup>*], *such that f*(*x*) = *<sup>∂</sup>x*(*q*1(*x*)*g*(*x*)), *and then it holds*

$$\mathcal{A}F\_i = \delta\_{\mathfrak{X}}(\varkappa \mathfrak{z}(\mathfrak{x})\delta\_{\mathfrak{X}}G)\_i + \mathcal{O}(h^4)\_i$$

*where Fi* = *f*(*xi*) *and Gi* = *g*(*xi*).

For any *ν*, *w* ∈ W*h*, we define the following inner products

$$
\langle \boldsymbol{\nu}, \boldsymbol{w} \rangle = \hbar \left[ \frac{1}{2} \nu\_0 w\_0 + \sum\_{i=1}^{M-1} \nu\_i w\_i + \frac{1}{2} \nu\_M w\_M \right], \quad \langle \boldsymbol{\nu}, \boldsymbol{w} \rangle\_1 = \hbar \sum\_{i=1}^{M-1} (\delta\_x w\_{i+1/2})(\delta\_x \boldsymbol{\nu}\_{i+1/2}),
$$

$$
\langle \boldsymbol{\nu}, \boldsymbol{w} \rangle\_{\varkappa\_2} = \hbar \sum\_{i=1}^{M-1} \varkappa \mathfrak{L}(\boldsymbol{x}\_{i+1/2})(\delta\_x w\_{i+1/2})(\delta\_x \boldsymbol{\nu}\_{i+1/2}), \quad \langle \boldsymbol{w}, \boldsymbol{\nu} \rangle\_{\mathcal{A}} = \langle \boldsymbol{w}, \boldsymbol{\nu} \rangle - \frac{\hbar^2}{12} \langle \boldsymbol{w}, \boldsymbol{\nu} \rangle.
$$

and their corresponding norms and semi-norms

$$\left| \left\langle w \right\rangle \right|^2 = \left\langle w, w \right\rangle\_{\prime} \quad \left| w \right|\_{1}^{2} = \left\langle w, w \right\rangle\_{1} \quad \left\| w \right\|\_{\mathcal{A}}^{2} = \left\langle w, w \right\rangle\_{\mathcal{A}} \quad \left| w \right|\_{1, \varkappa\_{2}}^{2} = \left\langle w, w \right\rangle\_{\varkappa\_{2}} \quad \left\| w \right\|\_{\infty} = \max\_{0 \le i \le M} |w\_{i}|.$$

Furthermore, assume that the coefficients satisfy

$$b\_0 \le q\_1(\mathbf{x}) \le b\_1, \quad b\_2 \le \varkappa\_2(\mathbf{x}) \le b\_3, \quad \left| \frac{q\_1'}{q\_1} \right| \le b\_3. \tag{13}$$

where all *bi* are positive constants.

**Lemma 2.6** ([25])**.** *For any grid function w* <sup>∈</sup> <sup>W</sup><sup>ˆ</sup> *h, it holds that*

$$\|w\| \le \frac{L}{\sqrt{6}} |w|\_{1\prime} \quad ||w||\_{\infty} \le \frac{\sqrt{L}}{2} |w|\_{1\prime} \quad \sqrt{\frac{2}{3}} \ ||w|| \le ||w||\_{\mathcal{A}} \le ||w||\\_{\mathcal{A}}$$

**Lemma 2.7** ([26])**.** *For any grid function w*1, *<sup>w</sup>*<sup>2</sup> <sup>∈</sup> <sup>W</sup><sup>ˆ</sup> *h, it holds that*

$$|\langle \mathcal{A}w, w \rangle| \ge \|w\|\_{\mathcal{A}}^2 - \frac{c\_3 h}{12} \left\| w \right\|^2, \quad |\langle \mathcal{A}w\_1, w\_2 \rangle| \le \frac{1}{2} \left[ \left\| w\_1 \right\|\_{\mathcal{A}}^2 + \left\| w\_2 \right\|\_{\mathcal{A}}^2 \right] + \frac{c\_3 h}{24} \left[ \left\| w\_1 \right\|^2 + \left\| w\_2 \right\|^2 \right].$$

**Lemma 2.8** ([8])**.** *For any <sup>h</sup>*(*t*) <sup>∈</sup> *<sup>C</sup>*3[0, *<sup>T</sup>*] *and <sup>γ</sup><sup>r</sup>* <sup>∈</sup> (0, 1], 0 <sup>≤</sup> *<sup>r</sup>* <sup>≤</sup> *m, such that <sup>γ</sup>*<sup>0</sup> <sup>&</sup>gt; *<sup>γ</sup>*<sup>1</sup> <sup>&</sup>gt; ··· <sup>&</sup>gt; *<sup>γ</sup>m*,*, the L*<sup>2</sup> − 1*<sup>σ</sup> formula has the following order of convergence*

*Mathematics* **2020**, *8*, 1696

$$\sum\_{r=0}^{m} p\_r \, {}\_0 D\_t^{\gamma\_r} h(t\_{k+r}) = \sum\_{n=0}^{k} \left( \sum\_{r=0}^{m} p\_r \frac{\tau^{-\gamma\_r}}{\Gamma(2-\gamma\_r)} \mathfrak{C}\_n^{(k+1,\gamma\_r)} \right) \left( h(t\_{k-n+1}) - h(t\_{k-n}) \right) + \mathcal{O}(\tau^{3-\gamma\_0}) \tag{14}$$

$$\hat{\mathbf{c}}\_{k} = \sum\_{n=0}^{k} \mathbb{C}\_{k-n}^{(k+1)} \left[ h(t\_{n+1}) - h(t\_n) \right] + \mathcal{O}(\tau^{3-\gamma\_0}).\tag{15}$$

**Lemma 2.9** ([24])**.** *For any h*(*t*) <sup>∈</sup> *<sup>C</sup>*3[0, *<sup>T</sup>*]*, we can obtain*

$$
\partial\_l h(t\_k) \cong \frac{1}{2\tau} \left[ (2\sigma + 1)h(t\_{k+1}) - 4\sigma h(t\_k) + (2\sigma - 1)h(t\_{k-1}) \right] \tag{16}
$$

$$k = \frac{\partial h}{\partial t}(t\_{k+\sigma}) + \mathcal{O}(\tau^2), \quad k \ge 1. \tag{17}$$

**Lemma 2.10** ([24])**.** *Suppose that* ·, ·∗ *is an inner product on* <sup>W</sup><sup>ˆ</sup> *<sup>h</sup> and* ·∗ *is a norm deduced by the inner product. For any grid functions w*0, *<sup>w</sup>*1,..., *wk*<sup>+</sup><sup>1</sup> <sup>∈</sup> <sup>W</sup><sup>ˆ</sup> *h, we have the following inequality*

$$\left\langle \sum\_{n=0}^{k} \hat{\mathfrak{L}}\_{k-n}^{(k+1)} \left[ w^{n+1} - w^{n} \right], w^{k+\upsilon} \right\rangle\_{\ast} \geq \frac{1}{2} \sum\_{n=0}^{k} \hat{\mathfrak{L}}\_{k-n}^{(k+1)} \left[ \left\| w^{n+1} \right\|\_{\ast}^{2} - \left\| w^{n} \right\|\_{\ast}^{2} \right],\tag{18}$$

$$\left< \partial\_l w^k, w^{k+\sigma} \right>\_\* \ge \frac{1}{4\pi} \left( \mathcal{E}^{k+1} - \mathcal{E}^k \right),\tag{19}$$

*where*

$$\mathcal{E}^{k+1} = (2\sigma + 1) \left\| w^{k+1} \right\|^2 - (2\sigma - 1) \left\| w^k \right\|^2 + (2\sigma^2 + \sigma - 1) \left\| w^{k+1} - w^k \right\|^2, \quad k \ge 0.$$

*Additionally, it holds that*

$$\mathcal{E}^{k+1} \ge \frac{1}{\sigma} \left\| w^{k+1} \right\|^2, \quad k \ge 0.$$

Let us initiate by order reduction for system (6) by letting *γ<sup>r</sup>* = *α<sup>r</sup>* − 1, (0 ≤ *r* ≤ *M*) and V(*x*, *t*) = *νt*(*x*, *t*), then

$$\frac{\partial^{\alpha\_r}\nu(\mathbf{x},t)}{\partial t^{\alpha\_r}} = \frac{\partial^{\gamma\_r}\mathcal{V}(\mathbf{x},t)}{\partial t^{\gamma\_r}}\nu$$

and so

$$
\partial\_t \left( \frac{\partial}{\partial \mathfrak{x}} \left( q\_1(\mathfrak{x}) \frac{\partial \nu}{\partial \mathfrak{x}} \right) \right) = \frac{\partial}{\partial \mathfrak{x}} \left( q\_1(\mathfrak{x}) \frac{\partial \mathcal{V}}{\partial \mathfrak{x}} \right) \dots
$$

Subsequently, the equivalent system to (6) after order reduction can be formulated as

$$\sum\_{\mathbf{r}=0}^{m} p\_{\mathbf{r}} \frac{\partial \gamma\_{\mathbf{r}} \mathcal{V}(\mathbf{x}, t)}{\partial t^{\gamma\_{\mathbf{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{20a}$$

$$
\partial\_t \left( \frac{\partial}{\partial x} \left( q\_1(\mathbf{x}) \frac{\partial \mathbf{v}}{\partial \mathbf{x}} \right) \right) = \frac{\partial}{\partial \mathbf{x}} \left( q\_1(\mathbf{x}) \frac{\partial \mathcal{V}}{\partial \mathbf{x}} \right), \tag{20b}
$$

with the following initial and boundary conditions

$$w(\mathbf{x},t) = \mathbf{\hat{r}}(\mathbf{x},t), \quad 0 \le \mathbf{x} \le L, \quad t \in [-s,0), \quad \mathcal{V}(\mathbf{x},0) = \boldsymbol{\hat{\psi}}(\mathbf{x}), \quad 0 \le \mathbf{x} \le L,\tag{20c}$$

$$v(0,t) = v(L,t) = 0, \quad 0 < t \le T,\tag{20d}$$

$$\mathcal{V}(0,t) = \mathcal{V}(L,t) = 0, \quad 0 < t \le T. \tag{20e}$$

#### *2.2. Compact Difference Scheme Construction*

Suppose that *<sup>ν</sup>t*(*x*, *<sup>t</sup>*) = <sup>V</sup>(*x*, *<sup>t</sup>*) <sup>∈</sup> **<sup>C</sup>**6,4 *<sup>x</sup>*,*<sup>t</sup>* ([0, *T*] × [0, *L*]), define the discretized functions

$$\mathbf{V}\_{i}^{k} = \mathcal{V}(\mathbf{x}\_{i\prime}t\_{k}), \quad \mathbf{V}\_{i}^{k} = \nu(\mathbf{x}\_{i\prime}t\_{k}), \; 0 \le i \le M, \; 0 \le k \le N.$$

Consider (20a) at (*xi*, *tk*<sup>+</sup>*σ*),, and then we get

$$\sum\_{r=0}^{m} p\_r \frac{\partial^{\gamma\_r} \mathcal{V}(\mathbf{x}\_i, t\_{k+\sigma})}{\partial t'^{\gamma\_r}} = \frac{\partial}{\partial \mathbf{x}} \left( q\_1(\mathbf{x}\_i) \frac{\partial \nu(\mathbf{x}\_i, t\_{k+\sigma})}{\partial \mathbf{x}} \right) + \tilde{f}\_i^{k+\sigma} \ 0 < i < M, \ 0 \le k \le N - 1,\tag{21}$$

such that

$$
\vec{f}\_i^k = \vec{f}(\nu(\mathbf{x}\_{i\prime}t\_k), \nu(\mathbf{x}\_{i\prime}t\_k - \mathbf{s}), \mathbf{x}\_{i\prime}t\_k).
$$

From Lemma 2.10, we conclude that

$$\sum\_{r=0}^{m} p\_r \frac{\partial^{\gamma\_r} \mathcal{V}(\mathbf{x}\_i, t\_{k+\tau})}{\partial t^{\gamma\_r}} = \sum\_{n=0}^{k} \hat{\mathbf{C}}\_{k-n}^{k+1} \left( \mathbf{V}\_i^{n+1} - \mathbf{V}\_i^n \right) + \mathcal{O}\left(\tau^{3-\gamma\_0}\right), \quad 0 < i < M, \ 0 \le k \le N-1,\tag{22}$$

and a direct expansion of Taylor type yields

$$
\tilde{f}\_i^{k+\sigma} = \tilde{F}\_i^{k+\sigma} + \mathcal{O}\left(\tau^2\right),
\tag{23}
$$

such that

$$\tilde{F}\_{i}^{k+\sigma} = \tilde{f}\left( (\sigma+1)\,\mathrm{V}\_{i}^{k} - \sigma \,\mathrm{V}\_{i}^{k-1}, \sigma \,\mathrm{V}\_{i}^{k+1-n\_{0}} + (1-\sigma)\,\mathrm{V}\_{i}^{k-n\_{0}}, \mathrm{x}\_{i}, t\_{k+\sigma} \right). \tag{24}$$

Acting the averaging operator A on both sides of (21), noticing Lemma 2.10 and using Taylor expansion, we arrive at

$$\sum\_{n=0}^{k} \mathcal{L}\_{k-n}^{k+1} \left( \mathcal{A} \mathbf{V}\_{l}^{n+1} - \mathcal{A} \mathbf{V}\_{l}^{n} \right) = \mathcal{A} \left[ \frac{\partial}{\partial \mathbf{x}} \left( q\_{1}(\mathbf{x}\_{l}) \frac{\partial v(\mathbf{x}\_{l}, t\_{k+\sigma})}{\partial \mathbf{x}} \right) \right] + \mathcal{A}\_{l}^{\mathbb{R}+\sigma} \ 0 < i < M, \ 0 \le k \le N - 1,\tag{25}$$

Next, by using Lemma 2.5 and (23), we obtain

$$\sum\_{n=0}^{k} \hat{\mathsf{C}}\_{k-n}^{k+1} \left( \mathcal{A} \mathbf{V}\_{i}^{n+1} - \mathcal{A} \mathbf{V}\_{i}^{n} \right) = \delta\_{\mathbf{x}} \left( \varkappa\_{2} \delta\_{\mathbf{x}} \mathbf{V} \right)\_{i}^{k+\sigma} + \mathcal{A} \vec{F}\_{i}^{k+\sigma} + \mathcal{S}\_{i}^{k+\sigma} \tag{26}$$

where a constant *c*˜0 exists in order that

$$|\mathcal{S}\_i^{k+\sigma}| \le \mathbb{E}\_0 \left(\tau^2 + h^4\right), \quad 1 \le i \le M - 1, \quad 0 \le k \le N - 1.$$

By considering (20b) at (*xi*, *t*1/2) and (*xi*, *tk*<sup>+</sup>*σ*), respectively, operating by A on both equations, we obtain by the aids of Taylor expansions and Lemmas 2.5 and 2.9, which

$$
\delta\_t \left( \delta\_{\mathbf{x}} \left( \varkappa\_2 \delta\_{\mathbf{x}} \mathbf{V} \right)\_i^{1/2} \right) = \delta\_{\mathbf{x}} \left( \varkappa\_2 \delta\_{\mathbf{x}} \mathbf{V} \right)\_i^{1/2} + s\_i^{1/2}, \quad 1 \le i \le M - 1,\tag{27}
$$

$$\partial\_{\hat{t}} \left( \delta\_{\mathbf{x}} \left( \varkappa\_2 \delta\_{\mathbf{x}} \mathbf{V} \right)\_{i}^{k+\sigma} \right) = \delta\_{\mathbf{x}} \left( \varkappa\_2 \delta\_{\mathbf{x}} \mathbf{V} \right)\_{i}^{k+\sigma} + s\_{i}^{k+\sigma}, \quad 1 \le i \le M-1, \quad 1 \le k \le N-1. \tag{28}$$

Moreover, there exists a constant *c*˜1 > 0, such that

$$|s\_i^{1/2}| \le \overline{c}\_1(\tau^2 + h^4), \quad 1 \le i \le M - 1,\tag{29}$$

$$|s\_i^{k+\sigma}| \le \varepsilon\_1 (\tau^2 + h^4), \quad 1 \le i \le M - 1, \quad 1 \le k \le N - 1. \tag{30}$$

By omitting the small terms in (26)–(28) and noticing the initial and boundary conditions, we construct a spatial fourth order difference scheme for problem (6), as follows

*Mathematics* **2020**, *8*, 1696

$$\sum\_{n=0}^{k} \mathfrak{C}\_{k-n}^{k+1} \left( \mathcal{A} \mathcal{V}\_{l}^{n+1} - \mathcal{A} \mathcal{V}\_{l}^{n} \right) = \delta\_{\mathbf{x}} \left( \varkappa\_{2} \delta\_{\mathbf{x}} \upsilon \right)\_{l}^{k+\sigma} + \mathcal{A} \mathcal{F}\_{l}^{k+\sigma}, \quad 1 \le i \le M-1, \quad 0 \le k \le N-1,\tag{31a}$$

$$
\delta\_l \left( \delta\_\mathbf{x} \left( \varkappa\_2 \delta\_\mathbf{x} \nu \right) \right)^{1/2} \right) = \delta\_\mathbf{x} \left( \varkappa\_2 \delta\_\mathbf{x} \mathcal{V} \right)^{1/2}, \quad 1 \le i \le M - 1,\tag{31b}
$$

$$\partial\_{\mathbb{I}}\left(\delta\_{\mathbb{X}}\left(\varkappa\_{2}\delta\_{\mathbb{X}}\nu\right)\_{l}^{k+\mathcal{O}}\right) = \delta\_{\mathbb{X}}\left(\varkappa\_{2}\delta\_{\mathbb{X}}\mathcal{V}\right)\_{l}^{k+\mathcal{O}}, \quad 1 \le i \le M-1, \quad 1 \le k \le N-1,\tag{31c}$$

$$\nu\_i^k = \mathbb{M}\_i^k, \quad 1 \le i \le M - 1, \quad -n\_0 \le k \le 0, \quad \mathcal{V}\_i^0 = \hat{\mathfrak{p}}(\mathbf{x}\_i), \quad 1 \le i \le M - 1,\tag{31d}$$

$$\nu\_0^k = \nu\_M^k = 0, \quad 0 \le k \le N,\tag{31e}$$

$$\mathcal{V}\_0^k = \mathcal{V}\_M^k = 0, \quad 0 \le k \le N. \tag{31}$$

where

$$\mathcal{F}\_{i}^{k+\sigma} = \tilde{f}\left( (\sigma+1)\,\nu\_{i}^{k} - \sigma\nu\_{i}^{k-1}, \sigma\nu\_{i}^{k+1-n\_{0}} + (1-\sigma)\nu\_{i}^{k-n\_{0}}, x\_{i}, t\_{k+\sigma} \right). \tag{32}$$
