**Case I** For arbitrary (*w*, *u*),(*w*∗, *u*∗) ∈ E and *x* ∈ [0, *x*1], we obtain

$$\begin{split} &|\, |\, \_\circ \mathcal{A}(w, u)(x) - \, \_\circ \mathcal{A}(w^\*, u^\*)(x)| \le |\rho(w(x)) - \rho(w^\*(x))| + \frac{1}{\Gamma(\varrho\_0)} \int\_0^x h\_0'(z)(h\_0(x) - h\_0(z))^{\varrho\_0 - 1} \\ &\times |\, \_\circ \mathcal{I}(z, u(z), w(z)) - \, \_\circ \mathcal{I}(z, u^\*(z), w^\*(z))| dz \\ &\le k\_\rho^\* |w(x) - w^\*(x)| + \frac{k\_f}{\Gamma(\varrho\_0)} \int\_0^x h\_0'(z)(h\_0(x) - h\_0(z))^{\varrho\_0 - 1} \\ &\times \left( |u - u^\*| + |w - w^\*| \right) dz \\ &\le \left( k\_\rho^\* + \frac{k\_f (h\_0(x\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \right) \|w - w^\*\| + \frac{k\_f (h\_0(x\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \|u - u^\*\|. \end{split} \tag{28}$$

Thus, we have

$$\begin{split} &|\mathcal{A}\_{1}(w,\boldsymbol{u})(\mathbf{x}) - \mathcal{A}\_{1}^{\*}(\boldsymbol{w}^{\*},\boldsymbol{u}^{\*})(\mathbf{x})| \\ &\leq \left(k\_{\rho}^{\*} + \frac{k\_{f}(h\_{0}(\mathbf{x}\_{1}) - h\_{0}(\mathbf{0}))^{\oplus\_{0}}}{\Gamma(\varrho\_{0} + 1)}\right) \|\boldsymbol{w} - \boldsymbol{w}^{\*}\| + \frac{k\_{f}(h\_{0}(\mathbf{x}\_{1}) - h\_{0}(\mathbf{0}))^{\oplus\_{0}}}{\Gamma(\varrho\_{0} + 1)} \|\boldsymbol{u} - \boldsymbol{u}^{\*}\|. \end{split} \tag{29}$$

Similarly

$$\begin{split} &|\mathcal{A}\_{2}(w,\boldsymbol{\mu})(\mathbf{x}) - \mathcal{A}\_{2}^{\prime}(w^{\*},\boldsymbol{\mu}^{\*})(\mathbf{x})| \\ &\leq \left(k\_{\boldsymbol{\phi}}^{\*} + \frac{k\_{\mathcal{F}}(h\_{0}(\mathbf{x}\_{1}) - h\_{0}(\mathbf{0}))^{\mathbf{e}\_{0}}}{\Gamma(\boldsymbol{\varrho}\_{0} + 1)}\right) ||\boldsymbol{\mu} - \boldsymbol{\mu}^{\*}|| + \frac{k\_{\mathcal{F}}(h\_{0}(\mathbf{x}\_{1}) - h\_{0}(\mathbf{0}))^{\mathbf{e}\_{0}}}{\Gamma(\boldsymbol{\varrho}\_{0} + 1)} ||\boldsymbol{w} - w^{\*}||. \end{split} \tag{30}$$

From (29) and (30), we have

$$\begin{split} & \| \mathcal{N} \left( w, u \right) - \mathcal{N} \left( w^\*, \mu^\* \right) \| \\ & \leq \left( k\_{\rho}^\* + k\_{\phi}^\* + 2 (k\_f + k\_{\mathcal{F}}) \frac{(h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \right) \left( \| w - w^\* \| + \| \mu - \mu^\* \| \right). \end{split} \tag{31}$$

**Case II**

For *x* ∈ (*x*k, *x*k+1], k = 1, . . . , ℵ, we have


Thus

$$\begin{split} \|\mathcal{N}\_{1}(\boldsymbol{w},\boldsymbol{u}) - \mathcal{N}\_{1}(\boldsymbol{w}^{\*},\boldsymbol{u}^{\*})\| &\leq k\_{f} \sum\_{i=0}^{\mathrm{k}} \frac{(h\_{i}(\boldsymbol{T}) - h\_{i}(\boldsymbol{x}\_{i}))^{\varrho\_{i}}}{\Gamma(\varrho\_{i} + 1)} \|\boldsymbol{u} - \boldsymbol{u}^{\*}\| \\ &+ \left(k\_{\rho}^{\*} + \mathbb{N}k\_{\mathcal{T}} + k\_{f} \sum\_{i=0}^{\mathrm{k}} \frac{(h\_{i}(\boldsymbol{T}) - h\_{i}(\boldsymbol{x}\_{i}))^{\varrho\_{i}}}{\Gamma(\varrho\_{i} + 1)}\right) \|\boldsymbol{w} - \boldsymbol{w}^{\*}\|. \end{split} \tag{33}$$

Similarly

$$\begin{split} \|\cdot\|\_{\mathcal{H}}(w,u) - \cdot\mathcal{H}(w^\*,u^\*)\| &\leq k\_{\mathcal{F}} \sum\_{i=0}^{\mathrm{k}} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \|w - w^\*\| \\ &+ \left(k\_{\Phi}^\* + \aleph k\_{\overline{\mathcal{T}}} + k\_{\mathcal{F}} \sum\_{i=0}^{\mathrm{k}} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)}\right) \|u - u^\*\|. \end{split} \tag{34}$$

From (33) and (34), we have

$$\left\|\left(\mathcal{N}\left(w,u\right)-\mathcal{N}\left(w^\*,u^\*\right)\right)\right\| \leq \left(k\_\rho^\* + k\_\phi^\* + \mathcal{N}\left(k\_\mathcal{T} + k\_\overline{\mathcal{T}}\right)\right) \tag{35}$$

$$\left(-2\left(k\_f + k\_{\mathcal{F}}\right)\sum\_{i=0}^\mathbf{k} \frac{\left(h\_i(T) - h\_i(\mathbf{x}\_i)\right)^{q\_i}}{\Gamma\left(q\_i + 1\right)}\right) \left(\left\|w - w^\*\right\| + \left\|u - u^\*\right\|\right).$$

Now if

$$\max(\chi\_{1\prime}\chi\_2) < 1,$$

where

$$\chi\_1 = k\_\rho^\* + k\_\phi^\* + 2(k\_f + k\_{\mathcal{F}}) \frac{(h\_0(\mathfrak{x}\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \lambda$$

and

$$\chi\_2 = k\_\rho^\* + k\_\phi^\* + \aleph(k\_\mathcal{T} + k\_{\overline{\mathcal{T}}}) + 2(k\_f + k\_{\mathcal{F}}) \sum\_{i=0}^\mathbf{k} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)}$$

then, N is strict contraction on E. It follows from Banach's contraction theorem that the impulsive FDE (1) has a unique solution on S.

#### **4. Stability Analysis of Problem** (1)

In this main section, we derive some results about stability analysis for the proposed problem (1). Prior to the proof of main results, we give definitions of Hyers–Ulam (H–U) stability and some remarks.

Consider an operator N : E→E, defined by

$$
\Box \mathcal{N}(w) = w; \quad w \in \mathcal{E}.\tag{36}
$$

**Definition 1.** *The solution w of problem* (36) *is H–U stable. If we find a constant* **C** > 0*, so that for any*  > 0 *and any solution w* ∈ E *of the inequality*

$$\left\{ \left| w - \mathcal{N}'(w) \right| \leq \epsilon \right. \tag{37}$$

*there exists unique solution w of Equation* (36) *in* E*, so that the following relation satisfies*

$$\|\overline{w} - w\| \le \mathsf{C}\epsilon.$$

**Definition 2.** *The solution of problem* (36) *is G–H–U stable if we find*

$$\theta : (0, \infty) \to (0, \infty), \theta(0) = 0$$

*so that for any solution of the inequality* (37)*, the following relation satisfies*

$$\|\overline{w} - w\| \le \mathbf{C}\theta(\epsilon).$$

**Remark 1.** *w is the solution in* E *for the inequality* (37)*, iff there exists a function* κ ∈ E *which is independent of solution* (*w*, *u*)*, so that for any t*


By Remark 1, we have the following perturbed problem

$$\begin{cases} ^cD\_{[x]}^{\varrho(x)}w(x) = f(x, u(x), w(x)) + \varkappa(x), \quad x \in \mathbb{S} = [0, T], x \neq x\_{i\prime} \\ w(0) = w\_0 + \rho(w), \\ \Delta w(\mathbf{x}\_i) = \mathbb{Z}\_i(w(\mathbf{x}\_i^{-})) + \varkappa\_{\mathbf{n}\_i}, \; i = 1, \ldots, m, \\ ^cD\_{[x]}^{\varrho(x)}u(\mathbf{x}) = \mathcal{F}(\mathbf{x}, w(\mathbf{x}), u(\mathbf{x})) + \varkappa(\mathbf{x}), \quad \mathbf{x} \in \mathbb{S} = [0, T], \mathbf{x} \neq \mathbf{x}\_{i\prime} \\ i = 1, \ldots, \aleph, \quad 0 < \varrho(\mathbf{x}) \le 1, \\ u(0) = u\_0 + \phi(u), \\ \Delta u(\mathbf{x}\_i) = \mathbb{Z}\_i(u(\mathbf{x}\_i^{-})) + \varkappa\_{\mathbf{n}\_i}, \; i = 1, \ldots, m. \end{cases} \tag{38}$$

**Lemma 2.** *The solution of the perturbed problem* (38) *satisfies the following relations*

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *w*(*x*) − *<sup>w</sup>*<sup>0</sup> <sup>+</sup> *<sup>ρ</sup>*(*w*) + <sup>1</sup> Γ(*-*0) *x* 0 *h* <sup>0</sup>(*z*)(*h*0(*x*) <sup>−</sup> *<sup>h</sup>*0(*z*))*-*<sup>0</sup>−<sup>1</sup> *<sup>f</sup>*(*z*, *<sup>u</sup>*(*z*), *<sup>w</sup>*(*z*))*dz* <sup>≤</sup> (*h*0(*x*1) <sup>−</sup> *<sup>h</sup>*0(0))*-*0 Γ(*-*<sup>0</sup> <sup>+</sup> <sup>1</sup>) , *if x* <sup>∈</sup> [0, *<sup>x</sup>*1], *. . . w*(*x*) − *w*<sup>0</sup> + *ρ*(*w*) + k ∑ *i*=1 I*iw*(*x*<sup>−</sup> *<sup>i</sup>* ) + k ∑ *i*=1 1 Γ(*<sup>i</sup>*−1) *xi xi*−<sup>1</sup> *h <sup>i</sup>*−1(*z*)(*hi*−1(*xi*) <sup>−</sup> *hi*−1(*z*))*<sup>i</sup>*−1−1 × *f*(*z*, *u*(*z*), *w*(*z*))*dz* + 1 Γ(k) *x x*k *h* <sup>k</sup>(*z*)(*h*k(*x*) <sup>−</sup> *<sup>h</sup>*k(*z*))<sup>k</sup>−<sup>1</sup> *<sup>f</sup>*(*z*, *<sup>u</sup>*(*z*), *<sup>w</sup>*(*z*))*dz* ≤ k + k ∑ *i*=0 (*hi*(*T*) <sup>−</sup> *hi*(*xi*))*i* Γ(*<sup>i</sup>* + 1) , *if x* ∈ (*x*k, *x*k+1], k = 1, . . . , ℵ, (39) *and* ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *u*(*x*) − *<sup>u</sup>*<sup>0</sup> <sup>+</sup> *<sup>φ</sup>*(*u*) + <sup>1</sup> Γ(*-*0) *x* 0 *h* <sup>0</sup>(*z*)(*h*0(*x*) <sup>−</sup> *<sup>h</sup>*0(*z*))*-*<sup>0</sup>−1F(*z*, *<sup>u</sup>*(*z*), *<sup>w</sup>*(*z*))*dz* <sup>≤</sup> (*h*0(*x*1) <sup>−</sup> *<sup>h</sup>*0(0))*-*0 Γ(*-*<sup>0</sup> <sup>+</sup> <sup>1</sup>) , *if x* <sup>∈</sup> [0, *<sup>x</sup>*1], *. . . u*(*x*) − *u*<sup>0</sup> + *φ*(*u*) + k ∑ *i*=1 I*iu*(*x*<sup>−</sup> *<sup>i</sup>* ) + k ∑ *i*=1 1 Γ(*<sup>i</sup>*−1) *xi xi*−<sup>1</sup> *h <sup>i</sup>*−1(*z*)(*hi*−1(*xi*) <sup>−</sup> *hi*−1(*z*))*<sup>i</sup>*−1−1 1 *x* (40)

$$\begin{split} & \mathbb{E} \left( \mathcal{F}(z, u(z), w(z)) dz + \frac{1}{\Gamma(\varrho\_{\mathbb{k}})} \int\_{\mathbf{x}\_{\mathbb{k}}}^{\mathbf{x}} h\_{\mathbf{k}}'(z) (h\_{\mathbf{k}}(\mathbf{x}) - h\_{\mathbf{k}}(\mathbf{z}))^{\varrho\_{\mathbb{k}} - 1} \mathcal{F}(z, u(z), w(z)) dz \right) \\ & \leq \left( \mathbb{k} + \sum\_{i=0}^{\mathbb{k}} \frac{(h\_{i}(T) - h\_{i}(\mathbf{x}\_{i}))^{\varrho\_{i}}}{\Gamma(\varrho\_{i} + 1)} \right) \varepsilon\_{i} \\ & \quad \text{if } \mathbf{x} \in (\mathbf{x}\_{\mathbb{k}}, \mathbf{x}\_{\mathbb{k}+1}], \mathbb{k} = 1, \dots, \mathbb{k}. \end{split}$$

**Proof.** The proof can be obtained by applying Lemma A2 repeatedly as in the proof of Lemma 1.

**Theorem 3.** *If* (*H*1), (*H*2) *and* (*H*6) *hold with the following condition*

$$\max(\chi\_1, \chi\_2) < 1,$$

*where*

$$\chi\_1 = k\_\rho^\* + k\_\phi^\* + 2(k\_f + k\_{\mathcal{F}}) \frac{(h\_0(\mathbf{x}\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \lambda$$

*and*

$$\chi\_2 = k\_\rho^\* + k\_\phi^\* + \aleph(k\_\mathcal{T} + k\_{\overline{\mathcal{T}}}) + 2(k\_f + k\_{\mathcal{F}}) \sum\_{i=0}^\mathbf{k} \frac{(h\_i(T) - h\_i(\chi\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)},$$

*then, problem* (1) *is H–U stable.*

**Proof.** Let *w*∗ be any solution of set of inequalities (37) and *w* be the unique solution of problem (1). Then, from integral Equations (8) and (39), we have

$$\begin{split} |w(\mathbf{x}) - w^\*(\mathbf{x})| &\leq |\rho(w(\mathbf{x})) - \rho(w^\*(\mathbf{x}))| + \frac{1}{\Gamma(\varrho\_0)} \int\_0^\mathbf{x} h\_0'(z) (h\_0(\mathbf{x}) - h\_0(\mathbf{z}))^{\varrho\_0 - 1} \\ &\geq |f(z, u(z), w(\mathbf{z})) - f(z, u^\*(z), w^\*(\mathbf{z}))| dz + \frac{1}{\Gamma(\varrho\_0)} \int\_0^\mathbf{x} h\_0'(z) (h\_0(\mathbf{x}) - h\_0(\mathbf{z}))^{\varrho\_0 - 1} |\varkappa(\mathbf{z})| dz \\ &\leq \left( k\_\rho^\* + \frac{k\_f (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \right) \|w - w^\*\| + \frac{k\_f (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \|u - u^\*\| \\ &+ \frac{\epsilon (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)}. \end{split} \tag{41}$$

Thus, for *x* ∈ [0, *x*1], we have

$$\begin{split} \|\|w - w^\*\|\| &\leq \left(k\_{\rho}^\* + \frac{k\_f(h\_0(\mathbf{x}\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)}\right) \|w - w^\*\| + \frac{k\_f(h\_0(\mathbf{x}\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \|u - u^\*\| \\ &+ \frac{\varepsilon(h\_0(\mathbf{x}\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)}. \end{split} \tag{42}$$

Similarly, for *x* ∈ [0, *x*1], we have

$$\begin{split} \|\|u - u^\*\|\| &\leq \left(k\_{\phi}^\* + \frac{k\_{\mathcal{F}}(h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)}\right) \|\|w - w^\*\|\| + \frac{k\_{\mathcal{F}}(h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \|\|u - u^\*\|\| \\ &+ \frac{\epsilon(h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)}. \end{split} \tag{43}$$

Adding (42) and (43), we have

$$\begin{split} & \| \boldsymbol{w} - \boldsymbol{w}^\* \| + \| \boldsymbol{u} - \boldsymbol{u}^\* \| \leq \left( k\_{\rho}^\* + \frac{k\_f (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \right) \| \boldsymbol{w} - \boldsymbol{u}^\* \| + \frac{k\_f (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \| \boldsymbol{u} - \boldsymbol{u}^\* \| \\ & + \frac{\epsilon (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \\ \leq \left( k\_{\rho}^\* + k\_{\phi}^\* + 2 (k\_f \\ & + k\_f \tau \frac{(h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \right) \left( \| \boldsymbol{w} - \boldsymbol{u} \| + \| \boldsymbol{w}^\* - \boldsymbol{u}^\* \| \right) + \frac{\epsilon (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)}. \end{split} \tag{44}$$

That implies

$$\begin{split} & \| \| w - w^\* \| + \| u - u^\* \| \le \left( k\_{\rho}^\* + k\_{\phi}^\* + 2(k\_f + k\_{\mathcal{F}}) \frac{(h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \right) \left( \| w - u \| + \| w^\* - u^\* \| \right) \\ & + \frac{\epsilon (h\_0(\mathbf{x}\_1) - h\_0(\mathbf{0}))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} . \end{split} \tag{45}$$

From which we obtain

$$\left\|(w,\mu)-(w^\*,\mu^\*)\right\| \le \frac{\frac{c(h\_0(x\_1)-h\_0(0))^{\epsilon\_0}}{\Gamma(\epsilon\_0+1)}}{1-\chi\_1},\tag{46}$$

where *χ*<sup>1</sup> = *k*∗ *<sup>ρ</sup>* + *k*<sup>∗</sup> *<sup>φ</sup>* <sup>+</sup> <sup>2</sup>(*<sup>k</sup> <sup>f</sup>* <sup>+</sup> *<sup>k</sup>*<sup>F</sup> )(*h*0(*x*1)−*h*0(0))*-*0 Γ(*-*<sup>0</sup>+1) is assumed to be less than one. By and large, for *x* ∈ (*x*k, *x*k+1], k = 1, . . . , ℵ, we have


Thus, we have

$$\begin{split} & \| |\boldsymbol{w} - \boldsymbol{w}^\* | \\ & \leq k\_f \sum\_{i=0}^{\mathbf{k}} \frac{\langle \boldsymbol{h}\_i(\boldsymbol{T}) - \boldsymbol{h}\_i(\mathbf{x}\_i) \rangle^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \| \boldsymbol{u} - \boldsymbol{u}^\* \| + \left( k\_\rho^\* + \aleph\_\mathcal{I} + k\_f \sum\_{i=0}^{\mathbf{k}} \frac{(\boldsymbol{h}\_i(\boldsymbol{T}) - \boldsymbol{h}\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \right) \| \boldsymbol{w} - \boldsymbol{w}^\* \| \\ & + \left( \mathbf{k} + \sum\_{i=0}^{\mathbf{k}} \frac{(\boldsymbol{h}\_i(\boldsymbol{T}) - \boldsymbol{h}\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \right) \varepsilon. \end{split} \tag{48}$$

Similarly, we have

$$\begin{split} & \| u - u^\* \| \\ & \le k\_{\mathcal{F}} \sum\_{i=0}^{\mathbf{k}} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \| w - w^\* \| + \left( k\_{\Phi}^\* + \aleph k\_{\mathcal{F}} + k\_{\mathcal{F}} \sum\_{i=0}^{\mathbf{k}} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \right) \| u - u^\* \| \\ & + \left( \mathbf{k} + \sum\_{i=0}^{\mathbf{k}} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \right) \varepsilon. \end{split} \tag{49}$$

From (48) and (49), we have

$$\begin{split} & \| \| w - w^\* \| + \| u - u^\* \| \le \left( k\_{\rho}^\* + k\_{\phi}^\* + \aleph(k\_{\mathcal{T}} + k\_{\overline{\mathcal{T}}} \big) + 2 \langle k\_f + k\_{\mathcal{F}} \rangle \sum\_{i=0}^{\mathrm{k}} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \right) \\ & \times \left( \| w - w^\* \| + \| u - u^\* \| \right) + \left( \mathrm{k} + \sum\_{i=0}^{\mathrm{k}} \frac{(h\_i(T) - h\_i(\mathbf{x}\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)} \right) \varepsilon. \end{split} \tag{50}$$

Which implies that

$$\left\| \left( w, u \right) - \left( w^\*, u^\* \right) \right\| \leq \frac{\left( \mathbb{k} + \sum\_{i=0}^{\mathbb{k}} \frac{\left( h\_i(T) - h\_i(\mathbf{x}\_i) \right)^{\rho\_i}}{\Gamma(\varrho\_i + 1)} \right) \epsilon}{1 - \chi\_2}. \tag{51}$$

Where *χ*<sup>2</sup> = *k*∗ *<sup>ρ</sup>* + *k*<sup>∗</sup> *<sup>φ</sup>* <sup>+</sup> <sup>ℵ</sup>(*k*<sup>I</sup> <sup>+</sup> *<sup>k</sup>*<sup>I</sup> ) + <sup>2</sup>(*<sup>k</sup> <sup>f</sup>* <sup>+</sup> *<sup>k</sup>*<sup>F</sup> ) <sup>∑</sup><sup>k</sup> *i*=0 (*hi*(*T*)−*hi*(*xi*))*i* Γ(*<sup>i</sup>*+1) is assumed to be less than one. Equivalently, (51) can be written as

$$\|\left(w,u\right)-\left(w^\*,u^\*\right)\|\leq \mathbf{C}\epsilon\_{\prime\prime}$$

where

$$\mathbf{C} = \frac{\left(\mathbf{k} + \sum\_{i=0}^{\mathbf{k}} \frac{(h\_i(T) - h\_i(x\_i))^{q\_i}}{\Gamma(q\_i + 1)}\right)}{1 - \left(k\_\rho^\* + k\_\phi^\* + \aleph(k\_\mathcal{T} + k\_\mathcal{T}) + 2(k\_f + k\_{\mathcal{F}})\sum\_{i=0}^{\mathbf{k}} \frac{(h\_i(T) - h\_i(x\_i))^{q\_i}}{\Gamma(q\_i + 1)}\right)}.$$

This shows that problem (1) is H–U stable.

**Lemma 3.** *By setting θ*() = **C**(); *θ*(0) = 0, *problem* (1) *becomes G–H–U stable.*

#### **5. Application and Discussion**

In this section, we apply our main results to the following numerical problem to verify the applications of the main results. We also plot graphs for its solution and functions *-* and *h* for illustration purposes.

#### **Example 1.**

$$\begin{cases} ^cD\_{[x]}^{\mathbf{e}(x)}w(\mathbf{x}) = \frac{e^{-\pi x}}{20} + \frac{(\mathbf{x} - \frac{1}{5})}{28} \left( |u(\mathbf{x})| + \sin(|w(\mathbf{x})|) \right), \\ \mathbf{x} \in [0, 1], \mathbf{x} \neq \mathbf{x}\_k, k = 1, 2, \dots, 9. \\ ^cD\_{[x]}^{\mathbf{e}(x)}u(\mathbf{x}) = \frac{e^{-\pi x}}{25} + \frac{(\mathbf{x} - \frac{1}{5})}{20} \left( |w(\mathbf{x})| + \cos(|u(\mathbf{x})|) \right), \\ \mathbf{x} \in [0, 1], \mathbf{x} \neq \mathbf{x}\_k, k = 1, 2, \dots, 9. \\ \Delta w \left( \mathbf{x}\_k \right) = \frac{1}{25} w(\mathbf{x}\_k^-), \quad \Delta u \left( \mathbf{x}\_k \right) = \frac{1}{40} u(\mathbf{x}\_k^-), \\ w(0) = \frac{w}{22 + |w|} + 0.025, \quad u(0) = \frac{u}{30 + |u|} + 1, \end{cases} \tag{52}$$

.

*where -* = <sup>1</sup> <sup>2</sup> , <sup>S</sup><sup>0</sup> = [0, <sup>1</sup> <sup>5</sup> ]*,* <sup>S</sup><sup>1</sup> = ( <sup>1</sup> <sup>5</sup> , 1]. *Set*

$$f(\mathbf{x}, u(\mathbf{x}), w(\mathbf{x})) = \frac{e^{-\pi \mathbf{x}}}{20} + \frac{(\mathbf{x} - \frac{1}{3})}{30} \left( |u(\mathbf{x})| + \sin(|w(\mathbf{x})|) \right); u, w \in \mathbb{R}^+,$$

*and*

$$\mathcal{F}(\mathbf{x}, u(\mathbf{x}), w(\mathbf{x})) = \frac{e^{-\pi x}}{25} + \frac{(\mathbf{x} - \frac{1}{5})}{20} \left( |w(\mathbf{x})| + \cos(|u(\mathbf{x})|) \right),$$

$$\mathcal{Z}\_i(w) = \frac{1}{50} w; w \in \mathbb{R}^+, i = 1, 2,$$

*and*

$$\rho(w) = \frac{w}{22 + |w|}, \quad \phi(u) = \frac{u}{30 + |u|}$$

*Assuming* ℵ = 2 *(*k = 1, 2*), we have*

$${}^{c}D\_{[x]}^{
\varrho(x)}w(x) = \begin{cases} {}^{c}D\_{[x]}^{
\varrho\_{0}J\_{0}}w(x), & 0 < x \le x\_{1}, \\ {}^{c}D\_{[x]}^{
\varrho\_{1}J\_{1}}^{
\varrho\_{1}}w(x), & x\_{1} < x \le x\_{2}, \\ {}^{c}D\_{[x]}^{
\varrho\_{2}J\_{2}}^{
\varrho\_{2}J\_{2}}w(x), & x\_{2} < x \le 1, \end{cases}$$

$$\, ^cD\_{[x]}^{\mathfrak{e}(x)}u(x) = \begin{cases} ^cD\_{[x]}^{\mathfrak{e}\_0J\_0}u(x), & 0 < x \le x\_{1\*} \\ ^cD\_{[x]}^{\mathfrak{e}\_1J\_1}u(x), & x\_1 < x \le x\_{2\*} \\ ^cD\_{[x]}^{\mathfrak{e}\_2J\_2}u(x), & x\_2 < x \le 1; \end{cases}$$

$$\varrho(x) = \begin{cases} \varrho\_0 = \frac{1}{4}, & 0 < x \le \frac{1}{3}, \\ \varrho\_1 = \frac{1}{3}, & \frac{1}{3} < x \le \frac{1}{2}, \\ \varrho\_2 = \frac{1}{2}, & \frac{1}{2} < x \le 1. \end{cases}$$

$$h(\mathbf{x}) = \begin{cases} h\_0(\mathbf{x}) = \frac{\mathbf{x}}{3}, & 0 < \mathbf{x} \le \frac{1}{3}, \\ h\_1(\mathbf{x}) = 2^{\mathbf{x}}, & \frac{1}{3} < \mathbf{x} \le \frac{1}{2}, \\ h\_2(\mathbf{x}) = e^{\mathbf{x}}, & \frac{1}{2} < \mathbf{x} \le 1. \end{cases}$$

*Let w*, *<sup>w</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> *and x* <sup>∈</sup> [0, 1]. *Then,*

$$\begin{split} & \left| f(\mathbf{x}, u(\mathbf{x}), w(\mathbf{x})) - f(\mathbf{x}, \overline{u}(\mathbf{x}), \overline{w}(\mathbf{x})) \right| \\ & \leq \frac{(\mathbf{x} - \frac{1}{5})}{28} \left( \left| |u(\mathbf{x}) - \overline{u}(\mathbf{x})| \right| + \left| \sin(|w(\mathbf{x})|) - \sin(\overline{w}(\mathbf{x})) \right| \right) \\ & \leq \frac{1}{35} \left( \left| |u(\mathbf{x}) - \overline{u}(\mathbf{x})| \right| + \left| \sin(|w(\mathbf{x})|) - \sin(\overline{w}(\mathbf{x})) \right| \right). \end{split} \tag{53}$$

*Similarly, we have*

$$\begin{split} & \left| \mathcal{F}(\mathbf{x}, u(\mathbf{x}), w(\mathbf{x})) - \mathcal{F}(\mathbf{x}, \overline{u}(\mathbf{x}), \overline{w}(\mathbf{x})) \right| \\ & \leq \frac{(\mathbf{x} - \frac{1}{5})}{20} \left( \left| |w(\mathbf{x}) - \overline{w}(\mathbf{x})| \right| + \left| \cos(|u(\mathbf{x})|) - \cos(\overline{u}(\mathbf{x})) \right| \right) \\ & \leq \frac{1}{25} \left( \left| |w(\mathbf{x}) - \overline{w}(\mathbf{x})| \right| + \left| \cos(|u(\mathbf{x})|) - \cos(\overline{u}(\mathbf{x})) \right| \right). \end{split} \tag{54}$$

*Using* (*H*1)*, from* (53) *and* (54)*, we obtain k <sup>f</sup>* = <sup>1</sup> <sup>35</sup> *and k*<sup>F</sup> <sup>=</sup> <sup>1</sup> <sup>25</sup> . *By* (*H*2)*,*

$$|\mathcal{Z}\_i(w) - \mathcal{Z}\_i(\overline{w})| \quad \le \quad \frac{1}{25}|w - \overline{w}|\_{\mathcal{A}}$$

$$|\mathbb{Z}\_i(\mu) - \mathbb{Z}\_i(\overline{\mu})| \le \frac{1}{40} |\mu - \overline{\mu}|.$$

*Using* (*H*2), *we get k*<sup>I</sup> <sup>=</sup> <sup>1</sup> <sup>25</sup> , *<sup>k</sup>*<sup>I</sup> <sup>=</sup> <sup>1</sup> 40 , *By* (*H*6), *we have*

$$\begin{aligned} |\rho(w) - \rho(\overline{w})| &= \left| \frac{w}{22 + |w|} - \frac{\overline{w}}{22 + |\overline{w}|} \right| \\ &\le \left| \frac{22|w - \overline{w}|}{(22 + |w|)(22 + |\overline{w}|)} \right| \le \frac{1}{22}|w - \overline{w}|, \\\\ |\phi(u) - \phi(\overline{u})| &= \left| \frac{u}{30 + |u|} - \frac{\overline{u}}{30 + |\overline{u}|} \right| \\ &\le \frac{30|u - \overline{u}|}{(30 + |u|)(30 + |\overline{u}|)} \le \frac{1}{30}|u - \overline{u}|. \end{aligned}$$

*Which implies k*∗ *<sup>φ</sup>* = <sup>1</sup> <sup>30</sup> . *Using the derived values, one may show that*

$$\max(\chi\_1, \chi\_2) < 1,$$

*where*

$$\chi\_1 = k\_\rho^\* + k\_\phi^\* + 2(k\_f + k\_{\mathcal{F}}) \frac{(h\_0(\varkappa\_1) - h\_0(0))^{\varrho\_0}}{\Gamma(\varrho\_0 + 1)} \prime$$

*and*

$$\chi\_2 = k\_\rho^\* + k\_\phi^\* + \aleph(k\_\mathcal{T} + k\_\mathcal{T}) + 2(k\_f + k\_{\mathcal{F}}) \sum\_{i=0}^\mathbb{k} \frac{(h\_i(T) - h\_i(\chi\_i))^{\varrho\_i}}{\Gamma(\varrho\_i + 1)}.$$

*Hence, by Theorem 2, the numerical problem* (52) *has a unique solution, and by Theorem 3, it is H–U stable. We have presented the piecewise graphs of function in Figure 1. The graph looks like a stair function. Moreover, the piecewise variable-order graphs for different pieces have been presented in Figure 2. The solution under the impulsive conditions and having piecewise variable-order has been plotted in Figure 3. The impulsive points are given as* 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9. *From the graph of solution, the crossover behaviors in the dynamics of the considered problem can be observed clearly at the given impulsive points. Hence, DEs with variable kernel have high flexibility due to the freedom of changing the kernel. This manuscript has a multiple stage structure. The problem investigated here has Caputo-type piecewise fractional-order derivative and a variable kernel. It can prove interesting for to readers and researchers working in this area.*

**Figure 1.** Plot for function in Example 1.

**Figure 2.** Plot for function *h* in Example 1.

**Figure 3.** Solution representation of problem (52) in Example 1.

### **6. Conclusions**

In this work, we have studied a coupled system of piecewise-order differential equations (DEs) with a variable kernel and impulsive conditions. The theoretical analysis is based on Scheafer's and Banach fixed-point theorems. For stability results, H–U's concept has been applied. The derived results have been applied to a numerical problem which illustrates the applicability of the main results. The contents of the paper generalize many results already studied in the literature. For the future, the reader should easily extend the results studied in [38,39] under the variable-order with a kernel of variable exponents. In addition, this concept can be extended to various problems of FDEs involving Caputo– Fabrizio or Atangana–Baleanu fractional differential operator with impulsive conditions and variable exponents.

**Author Contributions:** Conceptualization, A.A. (Arshad Ali); Methodology, K.J.A. and A.A. (Ahmad Aloqaily); Validation, H.A. and N.M.; Formal analysis, A.A. (Ahmad Aloqaily); Resources, K.J.A.; Data curation, H.A.; Writing—original draft, A.A. (Arshad Ali); Writing—review & editing, N.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors extend their appreciation to the Deanship of Scientific Research at King Khalid University for funding this work through large group Research Project under grant number RGP2/371/44. Ahmad Aloqaily, Nabil Mlaiki are thankful to Prince Sultan University for paying the APC and support through the TAS research lab.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

In this section, we give some definitions and preliminary results.

**Definition A1** ([6])**.** *The RL integral of fractional-order -*, *of function w*(*x*) *is given by*

$$\mathcal{Z}\_{a^{+}}^{\varrho}w(\mathbf{x}) = \frac{1}{\Gamma(\varrho)} \int\_{a}^{\mathbf{x}} (\mathbf{x} - \mathbf{s})^{\varrho - 1} w(\mathbf{s}) d\mathbf{s}.\tag{A1}$$

**Definition A2** ([24,40])**.** *The RL integral of fractional-order -*, *of function w*(*x*) *w.r.t h*(*x*) *is given by*

$$\mathcal{Z}\_{a^{+}}^{\varrho,h} w(\mathbf{x}) = \frac{1}{\Gamma(\varrho)} \int\_{a}^{\mathbf{x}} h'(\mathbf{s}) (h(\mathbf{x}) - h(\mathbf{s}))^{\varrho - 1} w(\mathbf{s}) d\mathbf{s};\tag{A2}$$

*the function h is increasing and differentiable such that h*(*x*) > 0, *for all x* > 0.

**Definition A3** ([6])**.** *The Caputo fractional derivative (CFD) of function w*(*x*) *is given by*

$$^cD\_{a^+}^\varrho w(\mathfrak{x}) = \mathcal{Z}\_{a^+}^{n-\varrho} w^{(n)}(\mathfrak{x}),\tag{A3}$$

*where n* − 1 < *-* < *n and w*(*n*)(*x*)=( *<sup>d</sup> dx* )*nw*(*x*).

**Definition A4** ([24,40])**.** *The CFD of function w*(*x*) *w.r.t h*(*x*) *is given by*

$${}^{c}D\_{a^{+}}^{\\\\\varrho,h}w(\mathbf{x}) = \mathcal{I}\_{a^{+}}^{n-\\\\\varrho,h}w\_{h}^{(n)}(\mathbf{x}),\tag{A4}$$

*where n* − 1 < *-* < *n and w*(*n*) *<sup>h</sup>* (*x*)=( <sup>1</sup> *h*(*x*) *d dx* )*nw*(*x*).

**Lemma A1** ([40])**.** *Let ϕ* ∈ *C*[*a*, *b*]*, a* < *b*, *so that the CFD exists. Then*

$$\prescript{c}{}{D}\_{a^{+}}^{\varrho,h} \prescript{}{a^{+}}{}{\varrho}(\mathfrak{x}) = \prescript{}{\varrho}(\mathfrak{x})\_{\prime}$$

*and*

$$\mathcal{T}\_{a^{+}}^{\varrho,h} \, ^{c}D\_{a^{+}}^{\varrho,h} \varphi(\mathfrak{x}) = \varrho(\mathfrak{x}) - \varrho(a)\_{\prime}$$

*for* 0 < *-* <sup>≤</sup> 1. *And cD-*,*h <sup>a</sup>*<sup>+</sup> *ϕ*(*x*) = 0 *if ϕ*(*x*) *is constant function.*

**Lemma A2** ([40])**.** *For -*∈ (0, 1], *the solution of the following problem*

$$\begin{aligned} \,^cD\_{a^+}^{\varrho,h} w(x) &=& \Phi(x),\\ w(a) &=& w\_0 \end{aligned} \tag{A5}$$

*is given by*

$$w(\mathbf{x}) = w\_0 + \frac{1}{\Gamma(\varrho)} \int\_a^\mathbf{x} h'(s) (h(\mathbf{x}) - h(s))^{\varrho - 1} \Phi(z) dz.$$

**Theorem A1.** *(Schaefer's fixed-point theorem) [41] Let* W *be a convex subset of a norm- linear space S with* 0 ∈ W *and let* B : W→W *is a completely continuous operator. Then the set* X = {*w* ∈ W : *w* = *ς*B*w*; 0 < *ς* < 1} *is either unbounded or* B *has a fixed point in* W*.*

#### **Appendix B**

The proof of Lemma 1 is received by using Lemma A2 for number of times. Assume *w* satisfies (5)–(7). If *x* ∈ [0, *x*1], then

$${}^{c}D\_{[x]}^{\\\otimes 0,h\otimes}w(\mathfrak{x}) = q(\mathfrak{x}), \\ [\mathfrak{x}] = 0.$$

Using Lemma A2, we get

$$w(x) = w\_0 + \rho(w) + \frac{1}{\Gamma(\varrho\_0)} \int\_0^x h\_0'(z) \left( h\_0(x) - h\_0(z) \right)^{\varrho\_0 - 1} \rho(z) dz.$$

This gives

$$w(\mathbf{x}\_1^{-}) = w\_0 + \rho(w) + \frac{1}{\Gamma(\varrho\_0)} \int\_0^{\mathbf{x}\_1} h\_0'(z) \left( h\_0(\mathbf{x}\_1) - h\_0(z) \right)^{\varrho\_0 - 1} \rho(z) dz.$$

Applying the impulse *w*(*x*− <sup>1</sup> ) = *<sup>w</sup>*(*x*<sup>+</sup> <sup>1</sup> ) − I1*w*(*x*<sup>−</sup> <sup>1</sup> ), we get

$$w(\mathbf{x}\_1^+) = w\_0 + \rho(w) + \frac{1}{\Gamma(\varrho\_0)} \int\_0^{\mathbf{x}\_1} h\_0'(z) (h\_0(\mathbf{x}\_1) - h\_0(z))^{\varrho\_0 - 1} \rho(z) dz + \mathcal{Z}\_1 w(\mathbf{x}\_1^-).$$

If *x* ∈ (*x*1, *x*2], then

$$\prescript{c}{}{D}\_{[x]}^{\varrho\_1, \dot{\eta}\_1} w(\mathfrak{x}) = \varrho(\mathfrak{x}), \,\, [\mathfrak{x}] = \mathfrak{x}\_1.$$

Using Lemma A2, we get

$$\begin{split} w(\mathbf{x}) &= \ &w(\mathbf{x}\_{1}^{+}) + \frac{1}{\Gamma(\varrho\_{1})} \int\_{x\_{1}}^{x} h\_{1}'(z) (h\_{1}(\mathbf{x}) - h\_{1}(z))^{\varrho\_{1} - 1} \varrho(z) dz \\ &= &w(\mathbf{x}\_{1}^{-}) + \mathcal{I}\_{1} w(\mathbf{x}\_{1}^{-}) + \frac{1}{\Gamma(\varrho\_{1})} \int\_{x\_{1}}^{x} h\_{1}'(z) (h\_{1}(\mathbf{x}) - h\_{1}(z))^{\varrho\_{1} - 1} \varrho(z) dz \\ &= &w\_{0} + \rho(w) + \frac{1}{\Gamma(\varrho\_{0})} \int\_{0}^{x\_{1}} h\_{0}'(z) (h\_{0}(\mathbf{x}\_{1}) - h\_{0}(z))^{\varrho\_{0} - 1} \varrho(z) dz \\ &+ \frac{1}{\Gamma(\varrho\_{1})} \int\_{x\_{1}}^{x} h\_{1}'(z) (h\_{1}(\mathbf{x}) - h\_{1}(z))^{\varrho\_{1} - 1} \varrho(z) dz + \mathcal{I}\_{1} w(\mathbf{x}\_{1}^{-}). \end{split}$$

This gives

$$\begin{aligned} w(\mathbf{x}\_2^-) &= \ &w\_0 + \rho(w) + \frac{1}{\Gamma(\varrho\_0)} \int\_0^{\mathbf{x}\_1} h\_0'(z) (h\_0(\mathbf{x}\_1) - h\_0(z))^{\varrho\_0 - 1} \varrho(z) dz \\ &+ \frac{1}{\Gamma(\varrho\_1)} \int\_{\mathbf{x}\_1}^{\mathbf{x}\_2} h\_1'(z) (h\_1(\mathbf{x}\_2) - h\_1(z))^{\varrho\_1 - 1} \varrho(z) dz + \mathcal{Z}\_1 w(\mathbf{x}\_1^-). \end{aligned}$$

Applying the impulse *w*(*x*− <sup>2</sup> ) = *<sup>w</sup>*(*x*<sup>+</sup> <sup>2</sup> ) − I2*w*(*x*<sup>−</sup> <sup>2</sup> ), we get

$$\begin{aligned} \left(w(\mathbf{x}\_2^+)\right)^\flat &= w\_0 + \rho(w) + \frac{1}{\Gamma(\varrho\_0)} \int\_0^{\mathbf{x}\_1} h\_0'(z) (h\_0(\mathbf{x}\_1) - h\_0(z))^{\varrho\_0 - 1} \varrho(z) dz \\ &+ \frac{1}{\Gamma(\varrho\_1)} \int\_{\mathbf{x}\_1}^{\mathbf{x}\_2} h\_1'(z) (h\_1(\mathbf{x}\_2) - h\_1(z))^{\varrho\_1 - 1} \varrho(z) dz + \mathcal{Z}\_1 w(\mathbf{x}\_1^-) + \mathcal{Z}\_2 w(\mathbf{x}\_2^-). \end{aligned}$$

If *x* ∈ (*x*2, *x*3], then

$$\, \, ^c D\_{[x]}^{\varrho\_2, \mathfrak{h}\_2} w(x) = \varrho(x), \, [x] = x\_2.$$

Using Lemma A2, we get

$$\begin{split} w(\mathbf{x}) &=& w(\mathbf{x}\_{2}^{+}) + \frac{1}{\Gamma(\varrho\_{2})} \int\_{x\_{2}}^{x} h\_{2}'(z) (h\_{2}(\mathbf{x}) - h\_{2}(z))^{\varrho\_{2}-1} \varrho(z) dz \\ &=& w(\mathbf{x}\_{2}^{-}) + \mathcal{Z}\_{2} w(\mathbf{x}\_{2}^{-}) + \frac{1}{\Gamma(\varrho\_{2})} \int\_{x\_{2}}^{x} h\_{2}'(z) (h\_{2}(\mathbf{x}) - h\_{2}(z))^{\varrho\_{2}-1} \varrho(z) dz \\ &=& w\_{0} + \rho(w) + \frac{1}{\Gamma(\varrho\_{0})} \int\_{0}^{x\_{1}} h\_{0}'(z) (h\_{0}(\mathbf{x}\_{1}) - h\_{0}(z))^{\varrho\_{0}-1} \varrho(z) dz \\ &+& \frac{1}{\Gamma(\varrho\_{1})} \int\_{x\_{1}}^{x\_{2}} h\_{1}'(z) (h\_{1}(\mathbf{x}\_{2}) - h\_{1}(z))^{\varrho\_{1}-1} \varrho(z) dz \\ &+& \frac{1}{\Gamma(\varrho\_{2})} \int\_{x\_{2}}^{x} h\_{2}'(z) (h\_{2}(\mathbf{x}) - h\_{2}(z))^{\varrho\_{2}-1} \varrho(z) dz + \mathcal{Z}\_{1} w(x\_{1}^{-}) + \mathcal{Z}\_{2} w(x\_{2}^{-}). \end{split}$$

This gives

$$\begin{split} w(\mathbf{x}\_{3}^{-}) &= \quad w\_{0} + \rho(w) + \frac{1}{\Gamma(\varrho\_{0})} \int\_{0}^{\mathbf{x}\_{1}} h\_{0}'(z) (h\_{0}(\mathbf{x}\_{1}) - h\_{0}(z))^{\varrho\_{0} - 1} \rho(z) dz \\ &+ \quad \frac{1}{\Gamma(\varrho\_{1})} \int\_{\mathbf{x}\_{1}}^{\mathbf{x}\_{2}} h\_{1}'(z) (h\_{1}(\mathbf{x}\_{2}) - h\_{1}(z))^{\varrho\_{1} - 1} \rho(z) dz \\ &+ \quad \frac{1}{\Gamma(\varrho\_{2})} \int\_{\mathbf{x}\_{2}}^{\mathbf{x}\_{3}} h\_{2}'(z) (h\_{2}(\mathbf{x}\_{3}) - h\_{2}(z))^{\varrho\_{2} - 1} \rho(z) dz \\ &+ \quad \mathcal{Z}\_{1} w(\mathbf{x}\_{1}^{-}) + \mathcal{Z}\_{2} w(\mathbf{x}\_{2}^{-}). \end{split}$$

Applying the impulse *w*(*x*− <sup>3</sup> ) = *<sup>w</sup>*(*x*<sup>+</sup> <sup>3</sup> ) − I3*w*(*x*<sup>−</sup> <sup>3</sup> ), we get

$$\begin{split} w(\mathbf{x}\_{3}^{+}) &= \; w\_{0} + \rho(w) + \frac{1}{\Gamma(\varrho\_{0})} \int\_{0}^{\mathbf{x}\_{1}} h\_{0}'(z) (h\_{0}(\mathbf{x}\_{1}) - h\_{0}(z))^{\varrho\_{0} - 1} \boldsymbol{\varrho}(z) dz \\ &+ \frac{1}{\Gamma(\varrho\_{1})} \int\_{\mathbf{x}\_{1}}^{\mathbf{x}\_{2}} h\_{1}'(z) (h\_{1}(\mathbf{x}\_{2}) - h\_{1}(z))^{\varrho\_{1} - 1} \boldsymbol{\varrho}(z) dz \\ &+ \frac{1}{\Gamma(\varrho\_{2})} \int\_{\mathbf{x}\_{2}}^{\mathbf{x}\_{3}} h\_{2}'(z) (h\_{2}(\mathbf{x}\_{3}) - h\_{2}(z))^{\varrho\_{2} - 1} \boldsymbol{\varrho}(z) dz \\ &+ \mathcal{Z}\_{1} w(\mathbf{x}\_{1}^{-}) + \mathcal{Z}\_{2} w(\mathbf{x}\_{2}^{-}) + \mathcal{Z}\_{3} w(\mathbf{x}\_{3}^{-}). \end{split}$$

Let

$$\begin{split} w(\mathbf{x}\_{\mathbf{k}}^{+}) &= \; w\_{0} + \mathcal{I}\_{1}w(\mathbf{x}\_{1}^{-}) + \mathcal{I}\_{2}w(\mathbf{x}\_{2}^{-}) + \mathcal{I}\_{3}w(\mathbf{x}\_{3}^{-}) + \dots + \mathcal{I}\_{k}w(\mathbf{x}\_{\mathbf{k}}^{-}) \\ &+ \int\_{0}^{\mathcal{T}} \frac{(T-z)^{\delta-1}}{\Gamma(\delta)} g(w(z)) dz + \frac{1}{\Gamma(\varrho\_{0})} \int\_{0}^{\mathbf{x}\_{1}} h\_{0}'(z) (h\_{0}(\mathbf{x}\_{1}) - h\_{0}(z))^{\varrho\_{0}-1} \boldsymbol{\varrho}(z) dz \\ &+ \quad \frac{1}{\Gamma(\varrho\_{1})} \int\_{\mathbf{x}\_{1}}^{\mathbf{x}\_{2}} h\_{1}'(z) (h\_{1}(\mathbf{x}\_{2}) - h\_{1}(z))^{\varrho\_{1}-1} \boldsymbol{\varrho}(z) dz + \frac{1}{\Gamma(\varrho\_{2})} \int\_{\mathbf{x}\_{2}}^{\mathbf{x}\_{3}} h\_{2}'(z) (h\_{2}(\mathbf{x}\_{3}) - h\_{2}(z))^{\varrho\_{2}-1} \boldsymbol{\varrho}(z) dz \\ &+ \quad \cdots + \frac{1}{\Gamma(\varrho\_{\mathbf{k}-1})} \int\_{\mathbf{x}\_{\mathbf{k}-1}}^{\mathbf{x}\_{\mathbf{k}}} h\_{\mathbf{k}-1}'(z) (h\_{\mathbf{k}-1}(\mathbf{x}\_{\mathbf{k}}) - h\_{\mathbf{k}-1}(z))^{\varrho\_{\mathbf{k}-1}-1} \boldsymbol{\varrho}(z) dz. \end{split}$$

Then, inductively, for *x* ∈ (*x*k, *x*k+1], we have

$$\prescript{\mathbf{c}}{}{D}\_{[\mathbf{x}]}^{\varrho\_{\mathbf{k}},\mu\_{\mathbf{k}}} w(\mathbf{x}) = \boldsymbol{\varrho}(\mathbf{x}), \, [\mathbf{x}] = \mathbf{x}\_{\mathbf{k}}.$$

Using Lemma A2, the solution becomes

$$\begin{split} w(\mathbf{x}) &= \quad w(\mathbf{x}\_{\mathbf{k}}^{+}) + \frac{1}{\Gamma(\varrho\_{\mathbf{k}})} \int\_{\mathbf{x}\_{\mathbf{k}}}^{\mathbf{x}} h\_{\mathbf{k}}'(z) (h\_{\mathbf{k}}(\mathbf{x}) - h\_{\mathbf{k}}(z))^{\varrho\_{\mathbf{k}} - 1} q(z) dz \\ &= \quad w\_{0} + \rho(w) + \sum\_{i=1}^{\mathbf{k}} \mathcal{Z}\_{i} w(\mathbf{x}\_{i}^{-}) \\ &+ \quad \sum\_{i=1}^{\mathbf{k}} \frac{1}{\Gamma(\varrho\_{i-1})} \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_{i}} h\_{i-1}'(z) (h\_{i-1}(\mathbf{x}\_{i}) - h\_{i-1}(z))^{\varrho\_{i-1} - 1} q(z) dz \\ &+ \quad \frac{1}{\Gamma(\varrho\_{\mathbf{k}})} \int\_{\mathbf{x}\_{\mathbf{k}}}^{\mathbf{x}} h\_{\mathbf{k}}'(z) (h\_{\mathbf{k}}(\mathbf{x}) - h\_{\mathbf{k}}(z))^{\varrho\_{\mathbf{k}} - 1} q(z) dz. \end{split}$$

Hence (4) holds. Conversely, let *w* satisfies the Equation (4). If *x* ∈ [0, *x*1], then *<sup>w</sup>*(0) = *<sup>w</sup>*0. Since *cD-*(*x*) [*x*] is the left inverse of I *-*(*x*) [*x*] thus using Lemma A1, we have

$$\prescript{\varepsilon}{}{D}\_{0}^{\varrho\_{0},h\_{0}}w(\mathfrak{x}) = \mathfrak{q}(\mathfrak{x}), \ \mathfrak{x} \in [0, \mathfrak{x}\_{1}].$$

If *<sup>x</sup>* <sup>∈</sup> [*x*k, *<sup>x</sup>*k+1), <sup>k</sup> <sup>=</sup> 1, ..., <sup>ℵ</sup>. Then for constant function *<sup>σ</sup>*(·), we have *cD-*(*x*) [*x*] *σ*(·) = 0. Thus

$${}^{c}D\_{[x]}^{\\\otimes \iota'h\_{\mathbb{K}}}w(\mathfrak{x}) = \mathfrak{q}(\mathfrak{x}), \text{ for each } \mathfrak{x} \in [\mathfrak{x}\_{\mathbb{K}}, \mathfrak{x}\_{\mathbb{k}+1}).$$

As well, we can simply infer that

$$w(\mathfrak{x}\_{\mathbb{k}}^{+}) - w(\mathfrak{x}\_{\mathbb{k}}^{-}) = \mathcal{Z}\_{\mathbb{k}}w(\mathfrak{x}\_{\mathbb{k}}^{-}), \mathbb{k} = 1, \dots, \aleph\_{\cdot}$$

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Fahad Alsidrani 1,2, Adem Kılıçman 2,3,***<sup>∗</sup>* **and Norazak Senu 2,3**

<sup>2</sup> Institute for Mathematical Research, Universiti Putra Malaysia, Serdang 43400, Selangor, Malaysia

<sup>3</sup> Department of Mathematics and Statistics, Faculty of Science, Universiti Putra Malaysia,


**Abstract:** In this research, three numerical methods, namely the variational iteration method, the Adomian decomposition method, and the homotopy analysis method are considered to achieve an approximate solution for a third-order time-fractional partial differential Equation (TFPDE). The equation is obtained from the classical (FW) equation by replacing the integer-order time derivative with the Caputo fractional derivative of order η = (0, 1] with variable coefficients. We consider homogeneous boundary conditions to find the approximate solutions for the bounded space variable *l* < χ < *L* and *l*, *L* ∈ R. To confirm the effectiveness of the proposed methods of non-integer order η, the computation of two test problems was presented. A comparison is made between the obtained results of the (VIM), (ADM), and (HAM) through tables and graphs. The numerical results demonstrate the effectiveness of the three numerical methods.

**Keywords:** fractional Fornberg–Whitham equation; approximate solution; partial differential equation; Riemann–Liouville derivatives; Caputo's derivatives; variational iteration method; Adomian decomposition method; homotopy analysis method

### **1. Introduction**

The concept of fractional partial differential equations (FPDEs) has been the focus of many studies and an essential topic in computational mathematics due to their various applications in scientific fields. The fractional derivative allows for a more accurate description of the diffusion process, considering the effects of long-range interactions and memory effects in most biological systems and phenomena in physics [1]. In recent years, researchers have demonstrated that many phenomena are successfully described by mathematical models of non-integer order using mathematical tools, for example, the Keller–Segel model for chemotaxis [1], fractional Riccati differential equations [2], and diffusion wave equations [3,4]. The cost of solving large nonlinear systems and related large linear systems after linearization can vary depending on various factors, including the complexity of the system and the method used for solving the system. It is worth noticing that different discretization methods for fractional diffusion equations (FDEs) have been proposed to solve a large linear system. In [5], Donatelli et al. have studied the diffusion equation, which arises in many applications that involves fractional derivatives in the case of variable coefficients (FDE). The proposed method was based on spectral analysis and structure preserving preconditioners for solving the (FDE). The method involves discretization in the space of the fractional diffusion equation, which leads to a linear system with coefficient matrices having a Toeplitz-like structure. In addition, they have shown that the variable coefficient matrix sequence belongs to the generalized locally Toeplitz (GLT) sequences. In [6], Donatelli et al. have implemented a finite volume (FV) method to discretize the space-fractional diffusion Equation (SFDE) with variable coefficients and obtain a large linear system resulting in a

**Citation:** Alsidrani, F.; Kılıçman, A.; Senu, N. Approximate Solutions for Time-Fractional Fornberg–Whitham Equation with Variable Coefficients. *Fractal Fract.* **2023**, *7*, 260. https:// doi.org/10.3390/fractalfract7030260

Academic Editors: Angelo B. Mingarelli, Leila Gholizadeh Zivlaei and Mohammad Dehghan

Received: 18 January 2023 Revised: 9 March 2023 Accepted: 11 March 2023 Published: 14 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Department of Mathematics, College of Science and Arts, Qassim University, Al Methnab 51931, Qassim, Saudi Arabia

sequence of coefficient matrices. The fractional derivative is considered in the Riesz-space fractional derivative. They showed that the resulting sequence of coefficient matrices belongs to the generalized locally Toeplitz (GLT) sequences. They also developed a good preconditioner and multigrid method to efficiently solve the obtained linear system of equations. In [7], Lin et al. developed fast algorithms for solving the linear systems that arise from the discretized time-dependent space-fractional diffusion Equation (SFDE) with non-constant coefficients. They also proved the convergence of two iteration schemes, one pre-smoother and the other post-smoother. In [8], Bu et al. presented a numerical method for solving a large linear system of the multi-term time-fractional advection–diffusion Equation (MTADE) using the finite element multigrid method. The method is based on the fractional derivative in the Riesz–Caputo sense, and the finite element approximation is considered in the space direction and time direction, respectively. They also discussed the stability and convergence of fully discrete schemes of (MTADE) in two situations.

The fractional derivative (FD) is even more significant in modeling real-life situations; for example, the fractional partial nonlinear Fornberg–Whitham (FPNFW) equation is a mathematical model that describes the evolution of nonlinear dispersive waves in fluid dynamics and the behavior of waves in plasmas.

Consider the nonlinear time-fractional Fornberg–Whitham equation

$$
\Psi^{\eta}\_{\zeta}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \psi\_{\chi\chi\zeta}(\boldsymbol{\chi},\boldsymbol{\zeta}) + \psi\_{\chi}(\boldsymbol{\chi},\boldsymbol{\zeta}) + \psi(\boldsymbol{\chi},\boldsymbol{\zeta})\psi\_{\chi}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \psi(\boldsymbol{\chi},\boldsymbol{\zeta})\psi\_{\chi\chi\chi}(\boldsymbol{\chi},\boldsymbol{\zeta}) - 3\psi\_{\chi}(\boldsymbol{\chi},\boldsymbol{\zeta})\psi\_{\chi\chi}(\boldsymbol{\chi},\boldsymbol{\zeta}) = 0 \tag{1}
$$

where ψ(χ, ζ) is the fluid velocity, 0 < η ≤ 1 is the order of fractional equation, ζ > 0 is the time, and χ is the spatial coordinate. In addition, when η = 1, Equation (1) is reduced to the original Fornberg–Whitham equation, which was first proposed by Whitham in 1967 for studying the qualitative behavior of wave breaking [9]. In 1978, Fornberg and Whitham [10]

obtained a peaked solution of the form <sup>ψ</sup>(χ, <sup>ζ</sup>) = <sup>K</sup>*exp* − 1/2|χ − 4ζ/3| , where K is an

arbitrary constant. In the literature, several mathematical methods have been implemented to obtain the approximate solutions of fractional differential equations, such as the Adomian decomposition method (ADM), variational iteration method (VIM), homotopy analysis method (HAM), homotopy perturbation method (HPM), Hermite wavelet method (HWM), optimal homotopy asymptotic method (OHAM), Shehu decomposition method (SDM), variational iteration transform method (VITM), Laplace decomposition method (LDM), direct power series method (DPSM), and others.

In [11], Kumar et al. solved the time-fractional Fornberg–Whitham equation involving the Atangana–Baleanu (AB) fractional derivative of non-integer order of the function ψ(χ, ζ) by using the Laplace decomposition method (LDM). This method is a mix of Adomian's decomposition method and the Laplace transform approach. The existence of the solution and the uniqueness of the solution of the nonlinear Fornberg–Whitham equation of fractional order were examined. In [12], Gupta and Singh used the homotopy perturbation method (HPM) to find the approximate numerical solution of the time-fractional Fornberg–Whitham equation where the derivatives are taken in the Caputo sense. In [13], Merdan et al. implemented a differential transformation method (DTM) to obtain an approximate analytical solution of the fractional Fornberg–Whitham equation. In [14], Alderremy et al. used the natural transform decomposition method (NTDM) to obtain the approximate numerical solution of the fractional Fornberg–Whitham equation in view of the Caputo operator. In [15], Fayçal and Omrani used two powerful techniques, namely the homotopy analysis method (HAM) and the Adomian's decomposition method (ADM) to obtain an approximate analytical solution of the nonlinear Fornberg–Whitham equation, where η = 1, and concluded that these methods have perfect accuracy and reductions in the size of calculations. In [16], Wang et al. combined He's (HPM) and the fractional complex transform to find an approximate solution to the nonlinear time-fractional Fornberg– Whitham equation. Recently, in [17], Sartanpara et al. used the p-Homotopy analysis Shehu transform method for the time-fractional Fornberg–Whitham equation with the derivative of the fractional-order in the Caputo sense to obtain the approximate analytical solution. In [18], Hijaz et al. numerically solved the Fornberg–Whitham classical type and modified type equations via the variational iteration algorithm-I. They used an auxiliary parameter to speed up the convergence rate to the exact solution. In [19], Shah et al. implemented modified techniques, namely the Shehu decomposition method (SDM) and the variational iteration transform method (VITM), to achieve an approximate analytical solution for the time-fractional Fornberg–Whitham equation. The fractional derivative is considered in the Caputo sense. In [20], Iqbal et al. successfully applied two modified methods to investigate the approximate solutions of the fractional Fornberg–Whitham equation. There was agreement between the numerical results obtained by the modified decomposition method (MDM) and modified variational iteration method (MVIM) involving fractional-order derivatives with Mittag–Leffler kernel.

In this paper, we consider the third-order time-fractional partial differential Equation (TFPDE) with variable coefficients

$$\begin{split} & \frac{\partial^{\mathsf{T}} \psi(\mathsf{X}, \mathsf{\zeta})}{\partial \mathsf{J}^{\mathsf{T}}} - \mathfrak{a}(\mathsf{x}) \frac{\partial^{3} \psi(\mathsf{X}, \mathsf{\zeta})}{\partial \mathsf{X}^{2} \partial \mathsf{U}} + 2\mathsf{e}(\mathsf{x}) \frac{\partial \psi(\mathsf{x}, \mathsf{\zeta})}{\partial \mathsf{X}} + \mathfrak{b}(\mathsf{x}) \psi(\mathsf{x}, \mathsf{\zeta}) \frac{\partial \psi(\mathsf{x}, \mathsf{\zeta})}{\partial \mathsf{X}} \\ & - \chi(\mathsf{x}) \psi(\mathsf{x}, \mathsf{\zeta}) \frac{\partial^{3} \psi(\mathsf{x}, \mathsf{\zeta})}{\partial \mathsf{X}^{3}} - \mathfrak{w}(\mathsf{x}) \frac{\partial \psi(\mathsf{x}, \mathsf{\zeta})}{\partial \mathsf{X}} \frac{\partial^{2} \psi(\mathsf{x}, \mathsf{\zeta})}{\partial \mathsf{X}^{2}} = 0, \qquad l < \mathsf{X} < L, \quad l, L \in \mathbb{R} \end{split} \tag{2}$$

with the initial and homogeneous boundary conditions

$$\begin{aligned} \Psi(\chi,0) &= \Phi(\chi) \\ \Psi(l,\zeta) &= 0, \quad \Psi(L,\zeta) = 0, \ \zeta &> 0 \end{aligned} \tag{3}$$

where α(*χ*), κ(*x*), β(*χ*), γ(*χ*), ω(*χ*) are the variable coefficients and *p* − <sup>1</sup> < η ≤ *p*,(*p* ∈ N) is a parameter describing the order of the time-fractional equation. The reason behind including variable coefficients in the time-fractional Fornberg–Whitham equation is that it becomes a more accurate model for the propagation of waves. This makes it a useful tool for studying a variety of phenomena, including fluid dynamics, plasma physics, and others. It is important to note that Equation (2) contains different interesting nonlinear equations. In particular, for η = 1, when β = 6 γ = −1 and α = κ = ω = 0 in Equation (2), we obtain the well-known Korteweg–de Vries Equation (KdV), as given in [21]

$$
\psi\_{\mathcal{L}} + 6\psi\psi\_{\mathcal{X}} + \psi\psi\_{\mathcal{XXX}} = 0, \quad \mathcal{X} \in \mathbb{R}, \quad \mathcal{L} > 0. \tag{4}
$$

When α = γ = 1, κ = 0, β = 4, and ω = 3 in Equation (2), we obtain the Degasperis– Procesi equation (DPE), as given in [22]

$$
\psi\_{\mathcal{L}} - \psi\_{\chi\chi\zeta} + 4\psi\psi\_{\chi} - \psi\psi\_{\chi\chi\chi} - 3\psi\_{\chi}\psi\_{\chi\chi} = 0, \quad \chi \in \mathbb{R}, \quad \zeta > 0. \tag{5}
$$

When α = ω = 1, β = 3, γ = 2 in Equation (2), and κ ∈ R is a parameter related to the critical shallow water speed, we obtain the Camassa–Holm Equation (CHE), as given in [21]

$$
\psi\_{\mathcal{L}} - \psi\_{\chi\chi\mathcal{L}} + 2\kappa\psi\_{\chi} + 3\psi\psi\_{\chi} - 2\psi\psi\_{\chi\chi\chi} - \psi\_{\chi}\psi\_{\chi\chi} = 0, \quad \chi \in \mathbb{R}, \quad \zeta > 0. \tag{6}
$$

When α = β = γ = 1, κ = <sup>1</sup> <sup>2</sup> and ω = 3 in Equation (2), we obtain the Fornberg–Whitham Equation (FWE), as given in [9]

$$
\psi\_{\mathcal{L}} - \psi\_{\chi\chi\mathcal{L}} + \psi\_{\chi} + \psi\psi\_{\chi} - \psi\psi\_{\chi\chi\chi} - \mathfrak{A}\psi\_{\chi}\psi\_{\chi\chi} = 0, \quad \chi \in \mathbb{R}, \quad \zeta > 0 \tag{7}
$$

We organize the paper as follows. Section 2 defines preliminary definitions and some properties of the Riemann–Liouville integral and Caputo fractional derivative. In Section 3, the analysis of (VIM) for the nonlinear fractional equation is established. In Section 4, the analysis of (ADM) for the nonlinear fractional equation is established. In Section 5, the analysis of (HAM) for the nonlinear fractional equation is established. Section 6

illustrates the methods for solving time-fractional partial differential equations (TFPDEs) with suitable initial conditions.

#### **2. Problem Formulation and Preliminaries**

In this section, we present the definitions of partial Riemann–Liouville integrals, partial Riemann–Liouville derivatives, and Caputo time-fractional derivatives with some properties of the Caputo fractional derivatives, which will be used later.

**Definition 1** ([23])**.** *Let* <sup>η</sup> <sup>∈</sup> (0, 1) *and* <sup>ψ</sup> <sup>∈</sup> *<sup>L</sup>*−1(*D*)*. The partial Riemann–Liouville fractional integrals of order* η *of a function* ψ(χ, ζ) *with respect to* ζ *are defined as*

$$\prescript{RL}{a}{\mathcal{T}}\_{\zeta}^{\eta}(\psi(\chi,\zeta)) = \frac{1}{\Gamma(\eta)} \int\_{a}^{\zeta} (\zeta - \xi)^{\eta - 1} \psi(\chi,\xi) \,\mathrm{d}\xi \tag{8}$$

*for almost all* (χ, ζ) ∈ *D and* Γ(η) *is the well-known Gamma function.*

**Definition 2** ([23])**.** *Let <sup>p</sup>* <sup>−</sup> <sup>1</sup> <sup>&</sup>lt; <sup>η</sup> <sup>≤</sup> *<sup>p</sup>*, *<sup>p</sup>* <sup>∈</sup> <sup>N</sup> *and* <sup>ψ</sup> <sup>∈</sup> *<sup>L</sup>*−1(*D*)*. The partial Riemann–Liouville fractional derivatives of order* η *of a function* ψ(χ, ζ) *with respect to* ζ *are defined as*

$$\, \_{\mathfrak{a}}\mathcal{D}\_{\zeta}^{\eta}(\psi(\chi,\zeta)) = \frac{\partial^{p}}{\partial \zeta^{p}} \, \_{\mathfrak{a}}\mathcal{D}\_{\zeta}^{p-\eta}(\psi(\chi,\zeta)) = \frac{\partial^{p}}{\partial \zeta^{p}} \left(\frac{1}{\Gamma(p-\eta)} \int\_{a}^{\zeta} (\zeta-\xi)^{p-\eta-1} \psi(\chi,\xi) d\xi\right) \tag{9}$$

*for almost all* (χ, ζ) ∈ *D.*

*For* <sup>ψ</sup> <sup>∈</sup> *<sup>C</sup>*μ*,* <sup>μ</sup> ≥ −1*,* <sup>γ</sup> <sup>&</sup>gt; <sup>−</sup><sup>1</sup> *and* <sup>η</sup>,<sup>β</sup> <sup>≥</sup> <sup>0</sup>*, the operator <sup>a</sup>*I<sup>η</sup> <sup>ζ</sup> *satisfies the following properties [24].*

*(1) <sup>a</sup>*I<sup>η</sup> <sup>ζ</sup> *<sup>a</sup>*I<sup>β</sup> <sup>ζ</sup> (ψ(ζ)) = *<sup>a</sup>*I<sup>β</sup> <sup>ζ</sup> *<sup>a</sup>*I<sup>η</sup> <sup>ζ</sup> (ψ(ζ)) = *<sup>a</sup>*Iη+<sup>β</sup> <sup>ζ</sup> (ψ(ζ)) *(2) <sup>a</sup>*I<sup>η</sup> <sup>ζ</sup> (<sup>ζ</sup> <sup>−</sup> *<sup>a</sup>*)<sup>γ</sup> <sup>=</sup> <sup>Γ</sup>(γ+1) <sup>Γ</sup>(η+γ+1)(<sup>ζ</sup> <sup>−</sup> *<sup>a</sup>*)η+<sup>γ</sup>

**Definition 3** ([25])**.** *Let p be the smallest integer that exceeds* η*, and the Caputo time-fractional derivative operator of order* η > 0 *of a function* ψ(χ, ζ) *is defined as*

$$\prescript{\mathsf{T}}{}{\mathcal{D}}\_{\zeta}^{\eta}(\mathsf{\boldsymbol{\upmu}}(\mathsf{y},\mathsf{\boldsymbol{\upzeta}})) = \begin{cases} \mathcal{T}\_{\zeta}^{p-\eta}\left(\frac{\partial^{p}\boldsymbol{\upphi}(\mathsf{y},\mathsf{\boldsymbol{\upzeta}})}{\partial\mathcal{J}^{p}}\right) = \frac{1}{\Gamma(p-\eta)}\int\_{0}^{\zeta}(\zeta-\xi)^{p-\eta-1}\frac{\partial^{p}\boldsymbol{\upphi}(\mathsf{y},\xi)}{\partial\mathcal{J}^{p}}d\xi, & p-1<\tau$$

*The operator* <sup>D</sup><sup>η</sup> <sup>ζ</sup> *satisfies the following properties [24]. Let* ζ > 0*, p* − 1 < η ≤ *p,* (*p* ∈ N)*, then*

$$\begin{aligned} (1) \quad & \mathcal{Z}\_{\boldsymbol{\zeta}}^{\boldsymbol{\eta}} \mathcal{D}\_{\boldsymbol{\zeta}}^{\boldsymbol{\eta}} (\boldsymbol{\psi}(\boldsymbol{\chi}, \boldsymbol{\zeta})) = \boldsymbol{\psi}(\boldsymbol{\chi}, \boldsymbol{\zeta}) - \sum\_{k=0}^{p-1} \frac{\zeta^{k}}{k!} \frac{\partial^{k}}{\partial \boldsymbol{\zeta}^{k}} \boldsymbol{\psi}(\boldsymbol{\chi}, \boldsymbol{\zeta}) \big|\_{\boldsymbol{\zeta} = 0} \\ (2) \quad & \mathcal{D}\_{\boldsymbol{\zeta}}^{\boldsymbol{\eta}} \mathcal{D}\_{\boldsymbol{\zeta}}^{\boldsymbol{\eta}} (\boldsymbol{\psi}(\boldsymbol{\chi}, \boldsymbol{\zeta})) = \boldsymbol{\psi}(\boldsymbol{\chi}, \boldsymbol{\zeta}) \end{aligned}$$

### **3. Analysis of Variational Iteration Method**

In this section, we discuss the (VIM) solution for the time-fractional partial differential Equation (TFPDE) with variable coefficients. This method can reduce the size of calculations and directly handle both linear and nonlinear equations, homogeneous or non-homogeneous [26,27].

We consider the following time-fractional partial differential equation

$$\mathcal{D}\_{\zeta}^{\eta} \psi(\chi, \zeta) + \mathcal{L} \psi(\chi, \zeta) + \mathcal{N} \psi(\chi, \zeta) = \mathcal{G}(\chi, \zeta) \tag{11}$$

where <sup>D</sup><sup>η</sup> <sup>ζ</sup> is the Caputo fractional derivative of order *p* − <sup>1</sup> < η ≤ *p*, (*p* ∈ N), L is the linear operator, N is the nonlinear operator, and G(χ, ζ) is a known analytical function. According to the variational iteration method [28,29], to solve the third-order time-fractional partial differential Equation (TFPDE) with variable coefficients in Equation (2), the correction functional can be constructed as follows

<sup>ψ</sup>*k*+1(χ, <sup>ζ</sup>) = <sup>ψ</sup>*k*(χ, <sup>ζ</sup>) + <sup>I</sup><sup>γ</sup> ζ λ Dη <sup>ζ</sup>ψ*k*(χ, <sup>ζ</sup>) + <sup>L</sup>ψ˜ *<sup>k</sup>*(χ, <sup>ζ</sup>) + <sup>N</sup> <sup>ψ</sup>˜ *<sup>k</sup>*(χ, <sup>ζ</sup>) <sup>=</sup> <sup>ψ</sup>*k*(χ, <sup>ζ</sup>) + <sup>1</sup> Γ(γ) ζ 0 (<sup>ζ</sup> <sup>−</sup> <sup>ξ</sup>)γ−1λ(ξ) *∂*ηψ*k*(χ, ξ) *<sup>∂</sup>*ξ<sup>η</sup> <sup>−</sup> <sup>α</sup>(χ) *∂*3ψ˜ *<sup>k</sup>*(χ, ξ) *<sup>∂</sup>*χ2*∂*<sup>ξ</sup> <sup>+</sup> <sup>2</sup>κ(χ) *∂*ψ˜ *<sup>k</sup>*(χ, ξ) *∂*χ + β(χ)ψ˜ *<sup>k</sup>*(χ, ξ) *∂*ψ˜ *<sup>k</sup>*(χ, ξ) *<sup>∂</sup>*<sup>χ</sup> <sup>−</sup> <sup>γ</sup>(χ)ψ˜ *<sup>k</sup>*(χ, <sup>ξ</sup>) *∂*3ψ˜ *<sup>k</sup>*(χ, ξ) *<sup>∂</sup>*χ<sup>3</sup> <sup>−</sup> <sup>ω</sup>(χ) *∂*ψ˜ *<sup>k</sup>*(χ, ξ) *∂*χ *∂*2ψ˜ *<sup>k</sup>*(χ, ξ) *∂*χ<sup>2</sup> dξ (12)

where <sup>I</sup><sup>γ</sup> <sup>ζ</sup> denotes the Riemann–Liouville integral operator of order γ = η + 1 − *p*, subject to the initial and boundary conditions in Equation (3), where λ is a general Lagrange multiplier which can be identified by variational theory, ψ*<sup>k</sup>* is the *k*th approximate solution and ψ˜ *<sup>k</sup>* is considered as a restricted variation, i.e., *δ*ψ˜ *<sup>k</sup>* = 0. This method first requires determining the Lagrange multiplier λ, and it can be easily identified as

$$\lambda(\xi) = \frac{(-1)^p (\xi - \zeta)^{p-1}}{(p-1)!} \tag{13}$$

where *p* is the highest order of the differential equation. By determining the value of Lagrange multiplier λ, the successive approximations ψ*k*+<sup>1</sup> will be calculated using the given initial function ψ0.

Making the above correction functional stationary and noticing that *δ*ψ˜ *<sup>k</sup>* = 0, we obtain

*<sup>δ</sup>*ψ*k*+1(χ, <sup>ζ</sup>) = *<sup>δ</sup>*ψ*k*(χ, <sup>ζ</sup>) + <sup>1</sup> Γ(γ) *δ* ζ 0 (<sup>ζ</sup> <sup>−</sup> <sup>ξ</sup>)γ−1λ(ξ) *∂*ηψ*k*(χ, ξ) *<sup>∂</sup>*ξ<sup>η</sup> <sup>−</sup> <sup>α</sup>(χ) *∂*3ψ˜ *<sup>k</sup>*(χ, ξ) *<sup>∂</sup>*χ2*∂*<sup>ξ</sup> <sup>+</sup> <sup>2</sup>κ(χ) *∂*ψ˜ *<sup>k</sup>*(χ, ξ) *∂*χ + β(χ)ψ˜ *<sup>k</sup>*(χ, ξ) *∂*ψ˜ *<sup>k</sup>*(χ, ξ) *<sup>∂</sup>*<sup>χ</sup> <sup>−</sup> <sup>γ</sup>(χ)ψ˜ *<sup>k</sup>*(χ, <sup>ξ</sup>) *∂*3ψ˜ *<sup>k</sup>*(χ, ξ) *<sup>∂</sup>*χ<sup>3</sup> <sup>−</sup> <sup>ω</sup>(χ) *∂*ψ˜ *<sup>k</sup>*(χ, ξ) *∂*χ *∂*2ψ˜ *<sup>k</sup>*(χ, ξ) *∂*χ<sup>2</sup> dξ <sup>=</sup> *<sup>δ</sup>*ψ*k*(χ, <sup>ζ</sup>) + <sup>1</sup> Γ(γ) *δ* ζ 0 (<sup>ζ</sup> <sup>−</sup> <sup>ξ</sup>)γ−1λ(ξ) *∂p*ψ*k*(χ, ξ) *∂*ξ*<sup>p</sup>* dξ (14)

This yields the Lagrange multipliers λ(ξ) = −1 for *p* = 1, and substituting this value of the Lagrange multiplier into the corrections functional Equation (12) gives the iteration formula for 0 < η ≤ 1

$$\begin{split} \Psi\_{k+1}(\mathbf{x},\boldsymbol{\zeta}) &= \Psi\_{k}(\mathbf{x},\boldsymbol{\zeta}) - \mathcal{Z}\_{\boldsymbol{\zeta}}^{\eta} \left( \frac{\partial^{\eta} \Psi\_{k}(\mathbf{x},\boldsymbol{\zeta})}{\partial \boldsymbol{\zeta}^{\eta}} - a(\boldsymbol{\chi}) \frac{\partial^{3} \Psi\_{k}(\mathbf{x},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{2} \partial \boldsymbol{\zeta}} + 2\kappa(\boldsymbol{\chi}) \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} \\ &+ \beta(\boldsymbol{\chi}) \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} - \gamma(\boldsymbol{\chi}) \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial^{3} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{3}} - a(\boldsymbol{\chi}) \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} \frac{\partial^{2} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{2}} \right) \end{split} (15)$$

Considering the given initial condition values, ψ0(χ, ζ) = φ(χ), and using this selection in Equation (15), we obtain the following successive approximations for *k* = 0, 1, ···

$$\begin{split} \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \Psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \mathcal{Z}\_{\zeta}^{\eta} \left( \frac{\partial^{\eta}\Psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial\boldsymbol{\zeta}^{\eta}} - a(\boldsymbol{\chi}) \frac{\partial^{3}\Psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}^{2}\partial\boldsymbol{\zeta}} + 2\kappa(\boldsymbol{\chi}) \frac{\partial\psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}} \\ &+ \beta(\boldsymbol{\chi})\Psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial\psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}} - \gamma(\boldsymbol{\chi})\psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial^{3}\psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}^{3}} - a(\boldsymbol{\chi}) \frac{\partial\psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}} \frac{\partial^{2}\psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}^{2}} \end{split} \tag{16}$$

$$\begin{split} \Psi\_{2}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \mathcal{T}\_{\boldsymbol{\zeta}}^{\eta} \left( \frac{\partial^{\eta} \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\zeta}^{\eta}} - \boldsymbol{\sigma}(\boldsymbol{\chi}) \frac{\partial^{3} \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{2} \partial \boldsymbol{\zeta}} + 2\kappa(\boldsymbol{\chi}) \frac{\partial \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} \\ &+ \beta(\boldsymbol{\chi}) \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} - \gamma(\boldsymbol{\chi}) \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial^{3} \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{3}} - \omega(\boldsymbol{\chi}) \frac{\partial \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} \frac{\partial^{2} \Psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{2}} \right) \tag{17} \\ &\vdots \end{split}$$

.

Thus, the corrections' functional Equation (15) will give a sequence of approximations <sup>ψ</sup>*k*(χ, <sup>ζ</sup>) = <sup>∑</sup>*k*−<sup>1</sup> *<sup>m</sup>*=<sup>0</sup> <sup>ψ</sup>*m*(χ, <sup>ζ</sup>).

Therefore, the solution of Equation (2) is given by

$$\psi(\chi,\zeta) = \lim\_{k \to \infty} \psi\_k(\chi,\zeta) \tag{18}$$

#### **4. Analysis of Adomian Decomposition Method**

In this section, we discuss the (ADM) solution for the time-fractional partial differential Equation (TFPDE) with variable coefficients. This method provides an analytical approximation to a rather wide class of nonlinear and stochastic equations without linearization, perturbation, closure approximations, or discretization methods resulting in massive numerical computation [30,31].

We consider the following time-fractional partial differential equation

$$
\mathcal{D}\_{\zeta}^{\eta} \psi(\chi, \zeta) + \mathcal{L} \psi(\chi, \zeta) + \mathcal{N} \psi(\chi, \zeta) = \mathcal{G}(\chi, \zeta) \tag{19}
$$

where <sup>D</sup><sup>η</sup> <sup>ζ</sup> is the Caputo fractional derivative of order *p* − 1 < η ≤ *p*, (*p* ∈ N), L is the linear operator, N is the nonlinear operator, and G(χ, ζ) is a known analytical function. To solve the third-order time-fractional partial differential Equation (TFPDE) with the variable coefficients shown in Equation (2), by Adomin decomposition method, we express this equation in the operator form as

$$\begin{cases} \mathcal{D}\_{\zeta}^{\eta} \psi(\chi, \zeta) - \mathfrak{a}(\chi) \mathcal{L}\_{\chi \chi \zeta} \psi(\chi, \zeta) + 2\mathfrak{a}(\chi) \mathcal{L}\_{\chi} \psi(\chi, \zeta) + \mathfrak{f}(\chi) \psi(\chi, \zeta) \mathcal{L}\_{\chi} \psi(\chi, \zeta) \\ - \chi(\chi) \psi(\chi, \zeta) \mathcal{L}\_{\chi \chi \chi} \psi(\chi, \zeta) - \mathfrak{a}(\chi) \mathcal{L}\_{\chi} \psi(\chi, \zeta) \mathcal{L}\_{\chi \chi} \psi(\chi, \zeta) = 0 \end{cases} \tag{20}$$

with the initial and boundary conditions shown in Equation (3), which are equivalent to

$$\mathcal{D}\_{\tilde{\zeta}}^{\eta} \psi(\chi, \zeta) - \alpha(\chi) \mathcal{L}\_{\chi \chi \zeta} \psi(\chi, \zeta) + 2\kappa(\chi) \mathcal{L}\_{\chi} \psi(\chi, \zeta) + \mathcal{N} \psi(\chi, \zeta) = 0 \tag{21}$$

Solving Equation (21) for <sup>D</sup><sup>η</sup> <sup>ζ</sup>ψ(χ, ζ), we obtain

$$\mathcal{D}\_{\zeta}^{\eta} \psi(\chi, \zeta) = \mathfrak{a}(\chi) \mathcal{L}\_{\chi \chi \zeta} \psi(\chi, \zeta) - 2\kappa(\chi) \mathcal{L}\_{\chi} \psi(\chi, \zeta) - \mathcal{N} \psi(\chi, \zeta) \tag{22}$$

where α(χ), κ(χ), β(χ), γ(χ), and ω(χ) are continuous functions, η is the parameter describing the order of the time-fractional derivative, the notations <sup>D</sup><sup>η</sup> <sup>ζ</sup> <sup>=</sup> *<sup>∂</sup>*<sup>η</sup> *<sup>∂</sup>*ζη , <sup>L</sup>χχζ <sup>=</sup> *<sup>∂</sup>*<sup>3</sup> *<sup>∂</sup>*χ2*∂*ζ, <sup>L</sup><sup>χ</sup> <sup>=</sup> *<sup>∂</sup> <sup>∂</sup>*<sup>χ</sup> , <sup>L</sup>χχ <sup>=</sup> *<sup>∂</sup>*<sup>2</sup> *<sup>∂</sup>*χ<sup>2</sup> , and <sup>L</sup>χχχ <sup>=</sup> *<sup>∂</sup>*<sup>3</sup> *<sup>∂</sup>*χ<sup>3</sup> are the symbolize of the linear operators, and N ψ(χ, ζ) = β(χ)ψ(χ, ζ)Lχψ(χ, ζ) −γ(χ)ψ(χ, ζ)Lχχχψ(χ, ζ) − ω(χ)Lχψ(χ, ζ)Lχχψ(χ, ζ) symbolizes the nonlinear operators. Applying the operator <sup>I</sup><sup>η</sup> <sup>ζ</sup> on both sides of Equation (22), with the basic properties of the operator <sup>D</sup><sup>η</sup> <sup>ζ</sup>, we obtain

$$\mathcal{Z}\_{\zeta}^{\eta} \mathcal{D}\_{\zeta}^{\eta} \psi(\chi, \zeta) = \mathcal{Z}\_{\zeta}^{\eta} \left[ \alpha(\chi) \mathcal{L}\_{\chi \chi \zeta} \psi(\chi, \zeta) - 2\kappa(\chi) \mathcal{L}\_{\chi} \psi(\chi, \zeta) - \mathcal{N} \psi(\chi, \zeta) \right] \tag{23}$$

$$\left|\psi(\mathbf{x},\boldsymbol{\zeta})-\sum\_{k=0}^{p-1}\frac{\boldsymbol{\zeta}^{k}}{k!}\frac{\partial^{k}}{\partial\boldsymbol{\zeta}^{k}}\psi(\mathbf{x},\boldsymbol{\zeta})\right|\_{\boldsymbol{\zeta}=0} = \mathcal{T}\_{\boldsymbol{\zeta}}^{\eta}\left[\boldsymbol{\alpha}(\boldsymbol{\chi})\mathcal{L}\_{\boldsymbol{\chi}\boldsymbol{\chi}\boldsymbol{\zeta}}\psi(\boldsymbol{\chi},\boldsymbol{\zeta})-2\kappa(\boldsymbol{\chi})\mathcal{L}\_{\boldsymbol{\chi}}\psi(\boldsymbol{\chi},\boldsymbol{\zeta})-\mathcal{N}\psi(\boldsymbol{\chi},\boldsymbol{\zeta})\right] \tag{24}$$

$$\Psi(\chi,\zeta) = \sum\_{k=0}^{p-1} \frac{\zeta^k}{k!} \frac{\partial^k}{\partial \zeta^k} \psi(\chi,\zeta) \bigg|\_{\zeta=0} + \mathcal{I}\_{\zeta}^{\eta} \left[ \mathfrak{a}(\chi) \mathcal{L}\_{\chi\chi\zeta} \psi(\chi,\zeta) - 2\kappa(\chi) \mathcal{L}\_{\chi} \psi(\chi,\zeta) - \mathcal{N} \psi(\chi,\zeta) \right] \tag{25}$$

Since the Adomain decomposition method is in the form of an infinite series

$$\Psi(\chi,\zeta) = \sum\_{m=0}^{\infty} \Psi\_m(\chi,\zeta) \tag{26}$$

and the nonlinear term N ψ(χ, ζ) can be decomposed into an infinite series of polynomials given by

$$\mathcal{N}\psi(\chi,\zeta) = \sum\_{m=0}^{\infty} \mathcal{A}\_m(\chi,\zeta) \tag{27}$$

where A*m*(χ, ζ) are the Adomain polynomials of ψ0, ψ1, ··· , ψ*<sup>m</sup>* defined by

$$\mathcal{A}\_m(\psi\_0, \psi\_1, \dots, \psi\_m) = \frac{1}{m!} \frac{d^m}{d\lambda^m} \left[ \mathcal{N} \left( \sum\_{j=0}^{\infty} \lambda^j \psi\_j \right) \right]\_{\lambda=0}, \quad m = 0, 1, \dots \tag{28}$$

Therefore, the first Adomain polynomials for N ψ(χ, ζ) are defined by

$$\begin{aligned} \mathcal{A}\_0 &= \mathcal{N}(\psi\_0) \\ \mathcal{A}\_1 &= \psi\_1 \left(\frac{d\mathcal{N}(\psi\_0)}{d\psi\_0}\right) \\ \mathcal{A}\_2 &= \psi\_2 \left(\frac{d\mathcal{N}(\psi\_0)}{d\psi\_0}\right) + \frac{\psi\_1^2}{2!} \left(\frac{d^2\mathcal{N}(\psi\_0)}{d\psi\_0^2}\right) \\ \mathcal{A}\_3 &= \psi\_3 \left(\frac{d\mathcal{N}(\psi\_0)}{d\psi\_0}\right) + \psi\_1 \psi\_2 \left(\frac{d^2\mathcal{N}(\psi\_0)}{d\psi\_0^2}\right) + \frac{\psi\_1^3}{3!} \left(\frac{d^3\mathcal{N}(\psi\_0)}{d\psi\_0^3}\right) \\ \mathcal{A}\_4 &= \psi\_4 \left(\frac{d\mathcal{N}(\psi\_0)}{d\psi\_0}\right) + \psi\_1 \psi\_3 \left(\frac{d^2\mathcal{N}(\psi\_0)}{d\psi\_0^2}\right) + \frac{\psi\_2^2}{2!} \left(\frac{d^2\mathcal{N}(\psi\_0)}{d\psi\_0^3}\right) + \frac{\psi\_1^3 \psi\_2}{2!} \left(\frac{d^3\mathcal{N}(\psi\_0)}{d\psi\_0^3}\right) + \frac{\psi\_1^3}{4!} N^{(4)}(\psi\_0) \\ \vdots \end{aligned} \tag{29}$$

Substituting Equations (26) and (27) into Equation (25), we obtain

$$\sum\_{m=0}^{\infty} \psi\_m(\boldsymbol{\chi}, \boldsymbol{\zeta}) = \sum\_{k=0}^{p-1} \frac{\boldsymbol{\zeta}^k}{k!} \frac{\partial^k}{\partial \boldsymbol{\zeta}^k} \psi(\boldsymbol{\chi}, \boldsymbol{\zeta}) \Big|\_{\boldsymbol{\zeta}=0} + \mathcal{Z}\_{\boldsymbol{\zeta}}^{\eta} \Big[ a(\boldsymbol{\chi}) \mathcal{L}\_{\boldsymbol{\chi}\boldsymbol{\zeta}\boldsymbol{\zeta}} \sum\_{m=0}^{\infty} \psi\_m(\boldsymbol{\chi}, \boldsymbol{\zeta}) - 2\kappa(\boldsymbol{\chi}) \mathcal{L}\_{\boldsymbol{\chi}} \sum\_{m=0}^{\infty} \psi\_m(\boldsymbol{\chi}, \boldsymbol{\zeta}) - \sum\_{m=0}^{\infty} \mathcal{A}\_m(\boldsymbol{\chi}, \boldsymbol{\zeta}) \Big] \tag{30}$$

The Adomain decomposition method transforms Equation (30) into a set of recursive relations given by

$$\begin{split} \psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \sum\_{k=0}^{p-1} \frac{\boldsymbol{\zeta}^{k}}{k!} \frac{\partial^{k}}{\partial \boldsymbol{\zeta}^{k}} \psi(\boldsymbol{\chi},\boldsymbol{\zeta}) \Big|\_{\boldsymbol{\zeta}=0} \\ \psi\_{m+1}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \mathcal{Z}\_{\boldsymbol{\zeta}}^{\text{T}} \Big[ \boldsymbol{\alpha}(\boldsymbol{\chi}) \mathcal{L}\_{\boldsymbol{\chi}\boldsymbol{\chi}\boldsymbol{\zeta}} \psi\_{\text{m}}(\boldsymbol{\chi},\boldsymbol{\zeta}) - 2\kappa(\boldsymbol{\chi}) \mathcal{L}\_{\boldsymbol{\chi}} \psi\_{\text{m}}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \mathcal{A}\_{\boldsymbol{m}}(\boldsymbol{\chi},\boldsymbol{\zeta}) \Big], \quad m \geq 0 \end{split} \tag{31}$$

Let the expression

$$\psi\_m(\chi,\zeta) = \sum\_{k=0}^{m-1} \psi\_k(\chi,\zeta) \tag{32}$$

be the *m*-term approximation of ψ. Using the above recursive relation Equation (31), we can obtain the first terms of (ADM) series solution for *m* = 0, 1, ···

$$\Psi\_1(\chi,\zeta) = \mathcal{Z}\_{\zeta}^{\eta} \left[ \alpha(\chi) \mathcal{L}\_{\chi\chi\zeta} \psi\_0(\chi,\zeta) - 2\kappa(\chi) \mathcal{L}\_{\chi} \psi\_0(\chi,\zeta) - \mathcal{A}\_0(\chi,\zeta) \right] \tag{33}$$

$$\Psi\_2(\chi,\zeta) = \mathcal{Z}\_{\zeta}^{\eta} \left[ \mathfrak{a}(\chi) \mathcal{L}\_{\chi\chi\zeta} \psi\_1(\chi,\zeta) - 2\kappa(\chi) \mathcal{L}\_{\chi} \psi\_1(\chi,\zeta) - \mathcal{A}\_1(\chi,\zeta) \right] \tag{34}$$

Therefore, the approximate solution is

. . .

$$\Psi(\chi,\zeta) = \lim\_{m \to \infty} \psi\_m(\chi,\zeta) \tag{35}$$

#### **5. Analysis of Homotopy Analysis Method**

In this section, we discuss the (HAM) solution for the time fractional partial differential Equation (TFPDE) with variable coefficients. Liao proposed a powerful and efficient method for nonlinear problems [32–34].

We consider the following time-fractional partial differential equation

$$
\mathcal{D}\_{\zeta}^{\eta} \psi(\chi, \zeta) + \mathcal{L} \psi(\chi, \zeta) + \mathcal{N} \psi(\chi, \zeta) = \mathcal{G}(\chi, \zeta) \tag{36}
$$

where <sup>D</sup><sup>η</sup> <sup>ζ</sup> is the Caputo fractional derivative of order *p* − 1 < η ≤ *p*, (*p* ∈ N), L is the linear operator, N is the nonlinear operator, and G(χ, ζ) is a known analytical function. To solve the third-order time-fractional partial differential Equation (TFPDE) with variable coefficients Equation (2) by homotopy analysis method, we consider the nonlinear operator

$$\begin{split} \mathcal{N}\left[\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})\right] &= \frac{\partial^{\mathsf{T}}\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})}{\partial\boldsymbol{\zeta}^{\mathsf{T}}} - a(\boldsymbol{\chi})\frac{\partial^{3}\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})}{\partial\boldsymbol{\chi}^{2}\partial\boldsymbol{\zeta}} + 2\kappa(\boldsymbol{\chi})\frac{\partial\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})}{\partial\boldsymbol{\chi}} + \beta(\boldsymbol{\chi})\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})\frac{\partial\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})}{\partial\boldsymbol{\chi}} \\ &- \gamma(\boldsymbol{\chi})\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})\frac{\partial^{3}\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})}{\partial\boldsymbol{\chi}^{3}} - \omega(\boldsymbol{\chi})\frac{\partial\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})}{\partial\boldsymbol{\chi}}\frac{\partial^{2}\boldsymbol{\varrho}(\boldsymbol{\chi},\boldsymbol{\zeta};\boldsymbol{\rho})}{\partial\boldsymbol{\chi}^{2}} = 0, \quad \boldsymbol{\zeta} > 0 \end{split} \tag{37}$$

and the linear operator

$$\mathcal{L}[\boldsymbol{\uprho}(\boldsymbol{\chi},\boldsymbol{\upzeta};\boldsymbol{\uprho})] = \mathcal{D}^{\eta}\_{\boldsymbol{\upzeta}}[\boldsymbol{\uprho}(\boldsymbol{\upchi},\boldsymbol{\upzeta};\boldsymbol{\uprho})] = \frac{\partial^{\eta}\,\boldsymbol{\uprho}(\boldsymbol{\upchi},\boldsymbol{\upzeta};\boldsymbol{\uprho})}{\partial\boldsymbol{\upzeta}^{\eta}}\tag{38}$$

subject to the initial and boundary conditions Equation (3), with the property <sup>D</sup><sup>η</sup> <sup>ζ</sup> (*k*) = 0, where *k* is the integration constant. According to Liao [33], we can construct the zero-order deformation equation

$$(1 - \rho) \mathcal{D}\_{\mathbb{Z}}^{\eta} \left[ \varphi(\chi, \mathbb{Q}; \rho) - \psi\_0(\chi, \mathbb{Q}) \right] = \rho \mathbb{M} \mathcal{H}(\chi, \mathbb{Q}) \mathcal{N} \left[ \varphi(\chi, \mathbb{Q}; \rho) \right] \tag{39}$$

where <sup>ρ</sup> <sup>∈</sup> [0, 1] is an embedding parameter, <sup>D</sup><sup>η</sup> <sup>ζ</sup> is an auxiliary linear operator, ϕ(χ, ζ; ρ) is a mapping function for ψ(χ, ζ), ψ0(χ, ζ) is an initial guess of ψ(χ, ζ), *h*¯ is a nonzero auxiliary parameter, and H(χ, ζ) is a nonzero auxiliary function. Obviously, for ρ = 0 and ρ = 1, we have

$$\begin{aligned} \varphi(\chi,\zeta;0) &= \psi\_0(\chi,\zeta) \\ \varphi(\chi,\zeta;1) &= \psi(\chi,\zeta) \end{aligned} \tag{40}$$

Thus, as ρ moves from 0 to 1, the solution ϕ(χ, ζ; ρ) varies from the initial guess ψ0(χ, ζ) to the solution ψ(χ, ζ). Expanding ϕ(χ, ζ; ρ) into the Taylor series with respect to the embedding parameter ρ, we obtain

$$\mathfrak{sp}(\chi,\zeta;\mathfrak{p}) = \mathfrak{\psi}\_0(\chi,\zeta) + \sum\_{m=1}^{\infty} \mathfrak{\psi}\_m(\chi,\zeta)\mathfrak{p}^m \tag{41}$$

where

$$\Psi\_{\rm{m}}(\chi,\zeta) = \frac{1}{m!} \frac{\partial^{m} \Phi(\chi,\zeta;\rho)}{\partial \rho^{m}} \Big|\_{\rho=0} \tag{42}$$

If the auxiliary linear operator, the initial guess, the auxiliary parameter *h*¯, and the auxiliary function are so properly chosen, the series Equation (41) converges at ρ = 1, then we have

$$\mathfrak{q}(\chi,\zeta;1) = \mathfrak{\psi}(\chi,\zeta) = \mathfrak{\psi}\_0(\chi,\zeta) + \sum\_{m=1}^{\infty} \mathfrak{\psi}\_m(\chi,\zeta) \tag{43}$$

which must be one of the solutions of the original nonlinear equation, as proven by Liao [35]. For ¯*h* = −1 and H(χ, ζ) = 1, Equation (39) becomes

$$\left[\mathfrak{l}(1-\mathfrak{p})\mathcal{D}\_{\zeta}^{\Pi}\right]\left[\mathfrak{q}(\chi,\zeta;\mathfrak{p})-\mathfrak{q}\mathfrak{q}(\chi,\zeta)\right]+\mathfrak{p}\mathcal{N}\left[\mathfrak{q}(\chi,\zeta;\mathfrak{p})\right]=0\tag{44}$$

According to Equation (42), the governing equation can be deduced from the zero-order deformation Equation (39).

Define the vector

$$\vec{\psi}\_m = \{ \psi\_0(\chi, \zeta), \psi\_1(\chi, \zeta), \dots, \psi\_m(\chi, \zeta) \}$$

By differentiating Equation (39) *m* number of times with respect to the embedding parameter ρ then setting ρ = 0, and finally dividing them by *m*!, we obtain the so-called *m*th-order deformation equation

$$\mathcal{D}\_{\vec{\zeta}}^{\eta}[\psi\_{m}(\chi,\zeta) - \chi\_{m}\psi\_{m-1}(\chi,\zeta)] = \hbar \mathcal{H}(\chi,\zeta)\mathcal{R}\_{m}[\vec{\psi}\_{m-1}(\chi,\zeta)] \tag{45}$$

where

$$\mathcal{R}\_{m}[\vec{\psi}\_{m-1}(\chi,\zeta)] = \frac{1}{(m-1)!} \frac{\partial^{m-1}\mathcal{N}[\varphi(\chi,\zeta;\rho)]}{\partial\rho^{m-1}}\Big|\_{\rho=0} \tag{46}$$

and

$$\chi\_m = \begin{cases} 0 & m \le 1 \\ 1 & m > 1 \end{cases} \tag{47}$$

Applying the operator <sup>I</sup><sup>η</sup> <sup>ζ</sup> on both sides of Equation (45), with the basic properties of the operator <sup>D</sup><sup>η</sup> <sup>ζ</sup>, we obtain the solution of the above *m*th-order deformation equation, with the assumption H(χ, ζ) = 1

$$
\Psi\_m(\chi,\zeta) = \chi\_m \Psi\_{m-1}(\chi,\zeta) + \hbar \mathcal{Z}\_{\zeta}^{\eta} \mathcal{R}\_m[\vec{\psi}\_{m-1}(\chi,\zeta)] \tag{48}
$$

where

.

$$\begin{split} \mathcal{R}\_{m}[\vec{\psi}\_{m-1}(\chi,\zeta)] &= \frac{\partial^{\eta}\psi\_{m-1}(\chi,\zeta)}{\partial\zeta^{\eta}} - \mathfrak{a}(\chi)\frac{\partial^{3}\psi\_{m-1}(\chi,\zeta)}{\partial\chi^{2}\partial\zeta} + 2\kappa(\chi)\frac{\partial\psi\_{m-1}(\chi,\zeta)}{\partial\chi} \\ &+ \sum\_{k=0}^{m-1} \left( \mathfrak{f}(\chi)\psi\_{k}(\chi,\zeta)\frac{\partial\psi\_{m-1-k}(\chi,\zeta)}{\partial\chi} - \mathfrak{v}(\chi)\psi\_{k}(\chi,\zeta)\frac{\partial^{3}\psi\_{m-1-k}(\chi,\zeta)}{\partial\chi^{3}} - \mathfrak{w}(\chi)\frac{\partial\psi\_{k}(\chi,\zeta)}{\partial\chi}\frac{\partial^{2}\psi\_{m-1-k}(\chi,\zeta)}{\partial\chi^{2}} \right) \end{split} \tag{49}$$

By using the above relation Equation (48) with the initial and boundary conditions Equation (3), we can obtain the first terms of the (HAM) series solution for *m* = 1, 2, ···

<sup>ψ</sup>1(χ, <sup>ζ</sup>) = *<sup>χ</sup>*1ψ0(χ, <sup>ζ</sup>) + *<sup>h</sup>*¯ <sup>I</sup><sup>η</sup> <sup>ζ</sup> <sup>R</sup>1[ ψ0(χ, ζ)] <sup>=</sup> *<sup>h</sup>*¯ <sup>I</sup><sup>η</sup> ζ *∂*ηψ0(χ, ζ) *<sup>∂</sup>*ζ<sup>η</sup> <sup>−</sup> <sup>α</sup>(χ) *∂*3ψ0(χ, ζ) *<sup>∂</sup>*χ2*∂*<sup>ζ</sup> <sup>+</sup> <sup>2</sup>κ(χ) *∂*ψ0(χ, ζ) *∂*χ + β(χ)ψ0(χ, ζ) *∂*ψ0(χ, ζ) *<sup>∂</sup>*<sup>χ</sup> <sup>−</sup> <sup>γ</sup>(χ)ψ0(χ, <sup>ζ</sup>) *∂*3ψ0(χ, ζ) *<sup>∂</sup>*χ<sup>3</sup> <sup>−</sup> <sup>ω</sup>(χ) *∂*ψ0(χ, ζ) *∂*χ *∂*2ψ0(χ, ζ) *∂*χ<sup>2</sup> (50) <sup>ψ</sup>2(χ, <sup>ζ</sup>) = *<sup>χ</sup>*2ψ1(χ, <sup>ζ</sup>) + *<sup>h</sup>*¯ <sup>I</sup><sup>η</sup> <sup>ζ</sup> <sup>R</sup>2[ ψ1(χ, ζ)] <sup>=</sup> <sup>ψ</sup>1(χ, <sup>ζ</sup>) + *<sup>h</sup>*¯ <sup>I</sup><sup>η</sup> ζ *∂*ηψ1(χ, ζ) *<sup>∂</sup>*ζ<sup>η</sup> <sup>−</sup> <sup>α</sup>(χ) *∂*3ψ1(χ, ζ) *<sup>∂</sup>*χ2*∂*<sup>ζ</sup> <sup>+</sup> <sup>2</sup>κ(χ) *∂*ψ1(χ, ζ) *∂*χ + β(χ)ψ0(χ, ζ) *∂*ψ1(χ, ζ) *<sup>∂</sup>*<sup>χ</sup> <sup>−</sup> <sup>γ</sup>(χ)ψ0(χ, <sup>ζ</sup>) *∂*3ψ1(χ, ζ) *<sup>∂</sup>*χ<sup>3</sup> <sup>−</sup> <sup>ω</sup>(χ) *∂*ψ0(χ, ζ) *∂*χ *∂*2ψ1(χ, ζ) *∂*χ<sup>2</sup> + β(χ)ψ1(χ, ζ) *∂*ψ0(χ, ζ) *<sup>∂</sup>*<sup>χ</sup> <sup>−</sup> <sup>γ</sup>(χ)ψ1(χ, <sup>ζ</sup>) *∂*3ψ0(χ, ζ) *<sup>∂</sup>*χ<sup>3</sup> <sup>−</sup> <sup>ω</sup>(χ) *∂*ψ1(χ, ζ) *∂*χ *∂*2ψ0(χ, ζ) *∂*χ<sup>2</sup> . . (51)

Therefore, we obtain an accurate approximation of Equation (2)

$$\Psi(\chi,\zeta) = \sum\_{k=0}^{m} \Psi\_k(\chi,\zeta) \tag{52}$$

#### **6. Applications and Results**

In this section, we apply (VIM), (ADM), and (HAM) to obtain the approximate solutions to the third-order time-fractional partial differential Equation (TFPDE) with variable coefficients and suitable initial conditions.

**Example 1.** *Consider the third-order time-fractional partial differential Equation (TFPDE) with variable coefficients*

$$\begin{split} & \frac{\partial^{\eta} \psi(\boldsymbol{\chi}, \boldsymbol{\zeta})}{\partial \boldsymbol{\zeta}^{\eta}} - \chi \frac{\partial^{3} \psi(\boldsymbol{\chi}, \boldsymbol{\zeta})}{\partial \chi^{2} \partial \boldsymbol{\zeta}} + 2\chi \frac{\partial \psi(\boldsymbol{\chi}, \boldsymbol{\zeta})}{\partial \chi} + \chi \psi(\boldsymbol{\chi}, \boldsymbol{\zeta}) \frac{\partial \psi(\boldsymbol{\chi}, \boldsymbol{\zeta})}{\partial \chi} \\ & - \chi \psi(\boldsymbol{\chi}, \boldsymbol{\zeta}) \frac{\partial^{3} \psi(\boldsymbol{\chi}, \boldsymbol{\zeta})}{\partial \chi^{3}} - \chi \frac{\partial \psi(\boldsymbol{\chi}, \boldsymbol{\zeta})}{\partial \chi} \frac{\partial^{2} \psi(\boldsymbol{\chi}, \boldsymbol{\zeta})}{\partial \chi^{2}} = 0, \quad \boldsymbol{\zeta} > 0, \quad 0 \le \boldsymbol{\chi} \le 1, \quad 0 < \eta \le 1 \end{split} \tag{53}$$

*and the initial and boundary conditions*

$$\begin{aligned} \Psi(\chi,0) &= \mathbf{e}^{\chi} \\ \Psi(0,\zeta) &= 0, \quad \Psi(1,\zeta) = 0, \quad \zeta > 0 \end{aligned} \tag{54}$$

Applying VIM: the iteration formula for Equation (53) can be constructed as

$$\begin{split} \Psi\_{k+1}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \Psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \mathcal{I}\_{\zeta}^{\eta} \left( \frac{\partial^{\eta} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\zeta}^{\eta}} - \chi \frac{\partial^{\boldsymbol{3}} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{2} \partial \boldsymbol{\zeta}} + 2\chi \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} \\ &+ \chi \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} - \chi \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial^{3} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{3}} - \chi \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}} \frac{\partial^{2} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\chi}^{2}} \end{split} \tag{55}$$

Considerin the given initial condition values, and using this selection in Equation (55), we obtain the following successive approximations

$$\begin{split} \psi\_{0}(\mathbf{x},\boldsymbol{\zeta}) &= \mathbf{e}^{\mathbf{x}} \\ \psi\_{1}(\mathbf{x},\boldsymbol{\zeta}) &= \mathbf{e}^{\mathbf{x}} (1 + \frac{\mathbf{x}(\mathbf{e}^{\mathbf{x}} - 2)\boldsymbol{\zeta}^{\eta}}{\Gamma(\eta + 1)}) \\ \psi\_{2}(\mathbf{x},\boldsymbol{\zeta}) &= \mathbf{e}^{\mathbf{x}} (1 + \frac{\mathbf{x}(\mathbf{e}^{\mathbf{x}} - 2)\boldsymbol{\zeta}^{\eta}}{\Gamma(\eta + 1)} + \frac{2\mathbf{x}(2\mathbf{q}\mathbf{x}^{\mathbf{x}} + 2\mathbf{e}^{\mathbf{x}} - \boldsymbol{\chi} - 2)\boldsymbol{\zeta}^{2\eta - 1}}{\Gamma(2\eta)} - \frac{4\mathbf{x}(-3\mathbf{q}\mathbf{e}^{2\mathbf{x}} - 4\mathbf{e}^{2\mathbf{x}} + 2\mathbf{q}\mathbf{e}^{\mathbf{x}} + 3\mathbf{e}^{\mathbf{x}} - \boldsymbol{\chi} - 1)\boldsymbol{\zeta}^{2\eta}}{\Gamma(2\eta + 1)} \\ &- \frac{\Gamma(2\eta + 1)\mathbf{x}\mathbf{e}^{\mathbf{x}}(2\mathbf{4}\mathbf{e}^{\mathbf{x}}\mathbf{y}^{\mathbf{x}} - 14\mathbf{q}^{2}\mathbf{e}^{2\mathbf{x}} + 52\mathbf{q}\mathbf{e}^{\mathbf{x}} - 23\mathbf{q}\mathbf{e}^{2\mathbf{x}} - 4\mathbf{q}^{2} + 12\mathbf{e}^{\mathbf{x}} - 4\mathbf{e}^{2\mathbf{x}} - 20\mathbf{y} - 8)\boldsymbol{\zeta}^{3\eta}}{\Gamma(3\eta + 1)\Gamma(\eta + 1)^{2}} \end{split} \tag{56}$$

Hence, the approximate solution for Equation (53) is

$$\begin{split} \psi \varphi(\mathbf{y},\boldsymbol{\zeta}) &= \mathbf{e}^{\chi} \left( 3 + \frac{2\mathbf{q}(\mathbf{e}^{\chi} - 2)\boldsymbol{\zeta}^{\eta}}{\Gamma(\eta + 1)} + \frac{2\chi(2\mathbf{q}\mathbf{e}^{\chi} + 2\mathbf{e}^{\chi} - \chi - 2)\boldsymbol{\zeta}^{2\eta - 1}}{\Gamma(2\eta)} - \frac{4\chi(-3\mathbf{q}\mathbf{e}^{2\chi} - 4\mathbf{e}^{2\chi} + 2\chi\mathbf{e}^{\chi} + 3\mathbf{e}^{\chi} - \chi - 1)\boldsymbol{\zeta}^{2\eta}}{\Gamma(2\eta + 1)} \\ &- \frac{\Gamma(2\eta + 1)\chi\mathbf{e}^{\chi}(24\mathbf{e}^{\chi}\mathbf{y}^{2} - 14\chi^{2}\mathbf{e}^{2\chi} + 52\chi\mathbf{e}^{\chi} - 23\chi\mathbf{e}^{2\chi} - 4\chi^{2} + 12\mathbf{e}^{\chi} - 4\mathbf{e}^{2\chi} - 20\chi - 8)\boldsymbol{\zeta}^{3\eta}}{\Gamma(3\eta + 1)\Gamma(\eta + 1)^{2}} \end{split} \tag{57}$$

Applying ADM: The recursive relations for Equation (53) can be constructed as

$$\begin{split} \Psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \sum\_{k=0}^{p-1} \frac{\boldsymbol{\zeta}^{k}}{k!} \frac{\partial^{k}}{\partial \boldsymbol{\zeta}^{k}} \boldsymbol{\Psi}(\boldsymbol{\chi},\boldsymbol{\zeta}) \Big|\_{\boldsymbol{\zeta}=0} \\ \boldsymbol{\Psi}\_{m+1}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \mathcal{Z}\_{\boldsymbol{\zeta}}^{\eta} \left[ \boldsymbol{\chi} \mathcal{L}\_{\chi\chi\zeta} \boldsymbol{\psi}\_{m}(\boldsymbol{\chi},\boldsymbol{\zeta}) - 2 \boldsymbol{\chi} \mathcal{L}\_{\chi} \boldsymbol{\psi}\_{m}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \mathcal{A}\_{m}(\boldsymbol{\chi},\boldsymbol{\zeta}) \right], \quad m \ge 0 \end{split} \tag{58}$$

Using the above recursive relations, we can obtain the first terms of the (ADM) series solution

$$\begin{split} \psi\_{0}(\mathbf{y},\boldsymbol{\zeta}) &= \mathbf{e}^{\mathbf{x}} \\ \psi\_{1}(\mathbf{y},\boldsymbol{\zeta}) &= \frac{\chi \mathbf{e}^{\mathbf{x}}(\mathbf{e}^{\mathbf{x}} - 2)\boldsymbol{\zeta}^{\eta}}{\Gamma(\eta + 1)} \\ \psi\_{2}(\mathbf{y},\boldsymbol{\zeta}) &= \chi \mathbf{e}^{\mathbf{x}} \left( \frac{(4\chi \mathbf{e}^{\mathbf{x}} + 4\mathbf{e}^{\mathbf{x}} - 2\chi - 4)\boldsymbol{\zeta}^{2\eta - 1}}{\Gamma(2\eta)} + \frac{(2\chi \mathbf{e}^{3\chi} + \mathbf{e}^{3\chi} - 4\chi \mathbf{e}^{2\chi} - 2\mathbf{e}^{2\chi} - 4\chi \mathbf{e}^{\chi} - 2\mathbf{e}^{\chi} + 4\chi + 4)\boldsymbol{\zeta}^{2\eta}}{\Gamma(2\eta + 1)} \right) \end{split} \tag{59}$$

Hence, the (ADM) series solution for Equation (53) is

$$\psi(\chi,\zeta) = \mathbf{e}^{\chi} \left( 1 + \frac{\chi(\mathbf{e}^{\chi} - 2)\zeta^{\eta}}{\Gamma(\eta + 1)} + \frac{\chi(4\chi\mathbf{e}^{\chi} + 4\mathbf{e}^{\chi} - 2\chi - 4)\zeta^{2\eta - 1}}{\Gamma(2\eta)} + \frac{\chi(2\chi\mathbf{e}^{3\chi} + \mathbf{e}^{3\chi} - 4\chi\mathbf{e}^{2\chi} - 2\mathbf{e}^{2\chi} - 4\chi\mathbf{e}^{\chi} - 2\mathbf{e}^{\chi} + 4\chi + 4)\zeta^{2\eta}}{\Gamma(2\eta + 1)} \right) \tag{60}$$

Applying HAM: the *m*th-order deformation equation for Equation (53) is given by

$$\begin{split} \Psi\_{m}(\mathbf{x},\boldsymbol{\zeta}) &= \chi\_{m}\Psi\_{m-1}(\mathbf{x},\boldsymbol{\zeta}) + \hbar\mathbb{Z}\_{\boldsymbol{\zeta}}^{\eta} \left( \frac{\partial^{\eta}\psi\_{m-1}(\mathbf{x},\boldsymbol{\zeta})}{\partial\boldsymbol{\zeta}^{\eta}} - \chi \frac{\partial^{3}\psi\_{m-1}(\mathbf{x},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}^{2}\partial\boldsymbol{\zeta}} + 2\chi \frac{\partial\psi\_{m-1}(\mathbf{x},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}} \\ &+ \sum\_{k=0}^{m-1} (\chi\psi\_{k}(\mathbf{y},\boldsymbol{\zeta})\frac{\partial\psi\_{m-1-k}(\mathbf{y},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}} - \chi\psi\_{k}(\mathbf{x},\boldsymbol{\zeta})\frac{\partial^{3}\psi\_{m-1-k}(\mathbf{y},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}^{3}} - \chi\frac{\partial\psi\_{k}(\mathbf{x},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}} \frac{\partial^{2}\psi\_{m-1-k}(\mathbf{y},\boldsymbol{\zeta})}{\partial\boldsymbol{\chi}^{2}} \end{split} \tag{61}$$

where

$$\chi\_{\mathfrak{M}} = \begin{cases} 0 & m \le 1 \\ 1 & m > 1 \end{cases} \tag{62}$$

Using the above relation Equation (61), we can obtain the first terms of (HAM) series solution

$$\begin{split} \psi \omega(\chi,\zeta) &= \mathbf{e}^{\chi} \\ \psi \omega\_{1}(\chi,\zeta) &= -\frac{\chi \mathbf{e}^{\chi} (\mathbf{e}^{\chi} - 2) \hbar \zeta^{\eta}}{\Gamma(\eta + 1)} \\ \psi \omega\_{2}(\chi,\zeta) &= \chi \mathbf{e}^{\chi} \left( -\frac{(\mathbf{e}^{\chi} - 2)(\hbar + 1) \hbar \zeta^{\eta}}{\Gamma(\eta + 1)} + \frac{4(3\chi \mathbf{e}^{2\chi} + 4\mathbf{e}^{2\chi} - 2\chi \mathbf{e}^{\chi} - 3\mathbf{e}^{\chi} + \chi + 1)\hbar^{2} \zeta^{2\eta}}{\Gamma(2\eta + 1)} + \frac{4((\chi + 1)\mathbf{e}^{\chi} - \frac{\chi}{2} - 1)\hbar^{2} \zeta^{2\tau - 1}}{\Gamma(2\eta)} \right) \end{split} \tag{63}$$

Hence, the (HAM) series solution for Equation (53) is

$$\psi(\chi,\zeta) = \mathbf{e}^{\chi} \left( 1 - \frac{\chi(\mathbf{e}^{\chi} - 2)\hbar^{2}\zeta^{\eta}}{\Gamma(\eta + 1)} + \frac{4\chi(3\chi\mathbf{e}^{2\chi} + 4\mathbf{e}^{2\chi} - 2\chi\mathbf{e}^{\chi} - 3\mathbf{e}^{\chi} + \chi + 1)\hbar^{2}\zeta^{2\eta}}{\Gamma(2\eta + 1)} + \frac{4\chi(\chi\mathbf{e}^{\chi} + \mathbf{e}^{\chi} - \frac{\chi}{2} - 1)\hbar^{2}\zeta^{2\tau - 1}}{\Gamma(2\eta)} \right) \tag{64}$$

**Example 2.** *Consider the third-order time-fractional partial differential Equation (TFPDE) with variable coefficients*

$$\begin{split} & \frac{\partial^{\eta} \psi(\mathbf{x}, \boldsymbol{\zeta})}{\partial \boldsymbol{\zeta}^{\eta}} - \chi^{2} \frac{\partial^{3} \psi(\mathbf{x}, \boldsymbol{\zeta})}{\partial \chi^{2} \partial \boldsymbol{\zeta}} + 2 \frac{\partial \psi(\mathbf{x}, \boldsymbol{\zeta})}{\partial \chi} + \psi(\mathbf{x}, \boldsymbol{\zeta}) \frac{\partial \psi(\mathbf{x}, \boldsymbol{\zeta})}{\partial \chi} \\ & - \psi(\mathbf{x}, \boldsymbol{\zeta}) \frac{\partial^{3} \psi(\mathbf{x}, \boldsymbol{\zeta})}{\partial \chi^{3}} - \frac{\partial \psi(\mathbf{x}, \boldsymbol{\zeta})}{\partial \chi} \frac{\partial^{2} \psi(\mathbf{x}, \boldsymbol{\zeta})}{\partial \chi^{2}} = 0, \quad \boldsymbol{\zeta} > 0, \quad 0 \le \boldsymbol{\chi} \le 1, \quad 0 < \boldsymbol{\eta} \le 1 \end{split} \tag{65}$$

*and the initial and boundary conditions*

$$\begin{aligned} \Psi(\chi,0) &= \chi^2\\ \Psi(0,\zeta) &= 0, \quad \Psi(1,\zeta) = 0, \quad \zeta > 0 \end{aligned} \tag{66}$$

Applying VIM: the iteration formula for Equation (65) can be constructed as

$$\begin{split} \Psi\_{k+1}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \Psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \mathcal{Z}\_{\boldsymbol{\zeta}}^{\eta} \left( \frac{\partial^{\eta} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \boldsymbol{\zeta}^{\eta}} - \chi^{2} \frac{\partial^{3} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \chi^{2} \partial \boldsymbol{\zeta}} + 2 \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \chi} \\ &+ \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \chi} - \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta}) \frac{\partial^{3} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \chi^{3}} - \frac{\partial \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \chi} \frac{\partial^{2} \psi\_{k}(\boldsymbol{\chi},\boldsymbol{\zeta})}{\partial \chi^{2}} \right) \end{split} \tag{67}$$

Considering the given initial condition values, and using this selection in Equation (66), we obtain the following successive approximations

$$\begin{split} \psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \boldsymbol{\chi}^{2} \\ \psi\_{1}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \boldsymbol{\chi}^{2} - \frac{(2\boldsymbol{\chi}^{3} - 8\boldsymbol{\chi})\boldsymbol{\zeta}^{\eta}}{\Gamma(\eta+1)} \\ \psi\_{2}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \boldsymbol{\chi}^{2} - \frac{(2\boldsymbol{\chi}^{3} - 8\boldsymbol{\chi})\boldsymbol{\zeta}^{\eta}}{\Gamma(\eta+1)} + \frac{12\boldsymbol{\chi}^{3}\boldsymbol{\zeta}^{2\eta-1}}{\Gamma(2\eta)} + \frac{(-10\boldsymbol{\chi}^{4} + 132\boldsymbol{\chi}^{2} - 32)\boldsymbol{\zeta}^{2\eta}}{\Gamma(2\eta+1)} + \frac{\Gamma(2\eta+1)(12\boldsymbol{\chi}^{5} - 304\boldsymbol{\chi}^{3} + 448\boldsymbol{\chi})\boldsymbol{\zeta}^{3\eta}}{\Gamma(3\eta+1)\Gamma(\eta+1)^{2}} \end{split} \tag{68}$$

Hence, the approximate solution for Equation (65) is

$$\psi(\chi\mathcal{L}) = 3\chi^2 - \frac{2(2\chi^3 - 8\chi)\zeta^\eta}{\Gamma(\eta + 1)} + \frac{12\chi^3\zeta^{2\eta - 1}}{\Gamma(2\eta)} + \frac{(-10\chi^4 + 132\chi^2 - 32)\zeta^{2\eta}}{\Gamma(2\eta + 1)} + \frac{\Gamma(2\eta + 1)(12\chi^5 - 304\chi^3 + 448\chi)\zeta^{2\eta}}{\Gamma(3\eta + 1)\Gamma(\eta + 1)^2} \tag{69}$$

Applying ADM: the recursive relations for Equation (65) can be constructed as

$$\begin{split} \Psi\_{0}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \sum\_{k=0}^{p-1} \frac{\boldsymbol{\zeta}^{k}}{k!} \frac{\partial^{k}}{\partial \boldsymbol{\zeta}^{k}} \Psi(\boldsymbol{\chi},\boldsymbol{\zeta}) \Big|\_{\boldsymbol{\zeta}=0} \\ \Psi\_{m+1}(\boldsymbol{\chi},\boldsymbol{\zeta}) &= \mathcal{Z}\_{\boldsymbol{\zeta}}^{\eta} \left[ \chi^{2} \mathcal{L}\_{\chi\chi\boldsymbol{\zeta}} \psi\_{m}(\boldsymbol{\chi},\boldsymbol{\zeta}) - 2 \mathcal{L}\_{\chi} \psi\_{m}(\boldsymbol{\chi},\boldsymbol{\zeta}) - \mathcal{A}\_{m}(\boldsymbol{\chi},\boldsymbol{\zeta}) \right], \quad m \ge 0 \end{split} \tag{70}$$

Using the above recursive relations, we can obtain the first terms of the (ADM) series solution

$$\begin{aligned} \psi\_0(\chi,\zeta) &= \chi^2\\ \psi\_1(\chi,\zeta) &= -\frac{(2\chi^3 - 8\chi)\zeta^\eta}{\Gamma(\eta + 1)}\\ \psi\_2(\chi,\zeta) &= \frac{4\left(3\chi^5 - 18\chi^3 + 3\chi^2 + 24\chi - 4\right)\zeta^{2\eta}}{\Gamma(2\eta + 1)} - \frac{12\chi^3\zeta^{2\eta - 1}}{\Gamma(2\eta)} \end{aligned} \tag{71}$$

Hence, the (ADM) series solution for Equation (65) is

$$\Psi(\chi,\zeta) = \chi^2 - \frac{(2\chi^3 - 8\chi)\zeta^\eta}{\Gamma(\eta + 1)} + \frac{4\left(3\chi^5 - 18\chi^3 + 3\chi^2 + 24\chi - 4\right)\zeta^{2\eta}}{\Gamma(2\eta + 1)} - \frac{12\chi^3\zeta^{2\eta - 1}}{\Gamma(2\eta)} \tag{72}$$

Applying HAM: the *m*th-order deformation equation for Equation (65) is given by

$$\begin{split} \Psi\_{\mathsf{m}}(\mathsf{y},\mathsf{\zeta}) &= \chi\_{\mathsf{m}}\psi\_{\mathsf{m}-1}(\mathsf{y},\mathsf{\zeta}) + \hbar\mathsf{I}\_{\mathsf{q}}^{\mathsf{m}} \Big( \frac{\partial^{\mathsf{m}}\psi\_{\mathsf{m}-1}(\mathsf{y},\mathsf{\zeta})}{\partial\mathsf{J}^{\mathsf{q}}} - \chi^{2} \frac{\partial^{\mathsf{3}}\psi\_{\mathsf{m}-1}(\mathsf{y},\mathsf{\zeta})}{\partial\mathsf{y}^{2}\partial\mathsf{J}} + 2\frac{\partial\psi\_{\mathsf{m}-1}(\mathsf{y},\mathsf{\zeta})}{\partial\mathsf{y}} \\ &+ \sum\_{k=0}^{\mathsf{m}-1} (\psi\_{\mathsf{k}}(\mathsf{y},\mathsf{\zeta})\frac{\partial\psi\_{\mathsf{m}-1-k}(\mathsf{y},\mathsf{\zeta})}{\partial\mathsf{y}} - \psi\_{\mathsf{k}}(\mathsf{y},\mathsf{\zeta})\frac{\partial^{\mathsf{3}}\psi\_{\mathsf{m}-1-k}(\mathsf{y},\mathsf{\zeta})}{\partial\mathsf{y}^{3}} - \frac{\partial\psi\_{\mathsf{k}}(\mathsf{y},\mathsf{\zeta})}{\partial\mathsf{y}} \frac{\partial^{2}\psi\_{\mathsf{m}-1-k}(\mathsf{y},\mathsf{\zeta})}{\partial\mathsf{y}^{2}} \Big) \end{split} \tag{73}$$

where

$$\chi\_m = \begin{cases} 0 & m \le 1 \\ 1 & m > 1 \end{cases} \tag{74}$$

Using the above relation Equation (73), we can obtain the first terms of (HAM) series solution

$$\begin{aligned} \Psi\_0(\chi,\zeta) &= \chi^2\\ \Psi\_1(\chi,\zeta) &= \frac{(2\chi^3 - 8\chi)\hbar\zeta^\eta}{\Gamma(\eta+1)}\\ \Psi\_2(\chi,\zeta) &= \frac{2\chi(\chi-2)(\chi+2)(\hbar+1)\hbar\zeta^\eta}{\Gamma(\eta+1)} + \frac{2(5\chi^4 - 66\chi^2 + 16)\hbar^2\zeta^{2\eta}}{\Gamma(2\eta+1)} - \frac{12\chi^3\hbar^2\zeta^{2\eta-1}}{\Gamma(2\eta)} \end{aligned} \tag{75}$$

Hence, the (HAM) series solution for Equation (65) is

$$\Psi(\mathbf{x},\boldsymbol{\zeta}) = \mathbf{x}^2 + \frac{2\chi(\chi-2)(\chi+2)(\hbar+2)\hbar\boldsymbol{\zeta}^{\eta}}{\Gamma(\eta+1)} + \frac{2(5\chi^4 - 6\epsilon\chi^2 + 16)\hbar^2\boldsymbol{\zeta}^{2\eta}}{\Gamma(2\eta+1)} - \frac{12\chi^3\hbar^2\boldsymbol{\zeta}^{2\eta-1}}{\Gamma(2\eta)}\tag{76}$$

#### **7. Conclusions**

This paper presents three numerical methods considered to achieve an approximate solution for a third-order time-fractional partial differential equation (TFPDE). The equation is obtained from the classical (FW) equation by replacing the integer-order time derivative with the Caputo fractional derivative of order η = (0, 1] with variable coefficients. The numerical results and graphs have been implemented using Maple 2022. In Figures 1–4 the graphical simulations for the approximate series solutions and the comparison in the form of absolute errors were presented to show the rate of change of the solutions when η = 0.75 and *h*¯ = −1. In Figures 5–8 we have shown the behavior of the solution with respect to the different values of η. As we can see in Tables 1 and 2, as ζ increases to 1, and χ increases to 1, the absolute errors slowly decrease. The agreement between the numerical results obtained by variational iteration method (VIM), Adomian decomposition method (ADM), and homotopy analysis method (HAM) involves fractional-order derivatives.

**Table 1.** This table shows the absolute errors for W = Abso.Error(VIM,ADM), U = Abso.Error(VIM,HAM), and V = Abso.Error(ADM,HAM) for Example 1.


**Table 1.** *Cont.* ζ χ *WUV* 00 2 2 0 0.1 0.1 2.219088227 2.152357150 0.066731077 0.2 0.2 2.628656823 2.243794656 0.384862167 0.3 0.3 3.627226957 2.269200185 1.358026772 0.4 0.4 6.084771319 2.229972343 3.854798976 η = 0.75 0.5 0.5 11.96807065 2.40088841 9.567182238 0.6 0.6 25.89804802 4.30750395 21.59054406 0.7 0.7 58.93797052 13.61183285 45.32613767 0.8 0.8 137.4703988 47.6301328 89.84026597 0.9 0.9 323.2302447 153.3892413 169.8410034 1 1 756.9909123 448.5531203 308.4377920

**Figure 1.** Graphical simulation of the second-level approximate solution ψ2(χ, ζ) when η = 0.75 and *h*¯ = −1 for Example 1.

**Figure 2.** Graphical simulation of the absolute error when η = 0.75 and ¯*h* = −1 for Example 1.

**Figure 6.** Plot of the absolute errors at χ = 0.1 for different values of η when ¯*h* = −1 for Example 1.

**Figure 7.** Plot of the first-level approximate solution ψ1(χ, ζ) at χ = 0.1 for different values of η when *h*¯ = −1 for Example 2.

**Table 2.** This table shows the absolute errors for W = Abso.Error(VIM,ADM), U = Abso.Error(VIM,HAM), and V = Abso.Error(ADM,HAM) for Example 2.


**Figure 8.** Plot of the absolute errors at χ = 0.1 for different values of η when ¯*h* = −1 for Example 2.

**Author Contributions:** Methodology, F.A.; Validation, A.K.; Investigation, A.K.; Writing—original draft, F.A.; Writing—review & editing, A.K.; Supervision, A.K. and N.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
