*Article* **Fractional Whitham–Broer–Kaup Equations within Modified Analytical Approaches**

#### **Rasool Shah1, Hassan Khan <sup>1</sup> and Dumitru Baleanu 2,3,\***


Received: 12 October 2019; Accepted: 5 November 2019; Published: 7 November 2019

**Abstract:** The fractional traveling wave solution of important Whitham–Broer–Kaup equations was investigated by using the q-homotopy analysis transform method and natural decomposition method. The Caputo definition of fractional derivatives is used to describe the fractional operator. The obtained results, using the suggested methods are compared with each other as well as with the exact results of the problems. The comparison shows the best agreement of solutions with each other and with the exact solution as well. Moreover, the proposed methods are found to be accurate, effective, and straightforward while dealing with the fractional-order system of partial differential equations and therefore can be generalized to other fractional order complex problems from engineering and science.

**Keywords:** q-Homotopy analysis transform method; Natural decomposition method; Whitham–Broer–Kaup equations; Caputo derivative

#### **1. Introduction**

The modern, broadly considered concept of fractional calculus was developed from a question raised by L'Hospital to Gottfried Wilhelm Leibniz in 1695. L'Hospital insisted on knowing about the outcome of the derivative of order *α* = <sup>1</sup> <sup>2</sup> , which laid down the foundation of a powerful fractional calculus [1,2]. Since then, the new theory of fractional calculus has gained the full attention of mathematicians, physicists, biologists, engineers, and economists in many areas of applied science. In modern decades, researchers have recognized that fractional-order differential equations contributed, in a natural way, to the study of different physical problems, such as diffusion processes, signal processing, viscoelastic systems, control processing, fractional stochastic systems, biology and ecology, quantum mechanics, wave theory, biophysics, and other research fields [3,4].

Partial differential equations (PDEs) involving non-linearities explain different phenomena in applied sciences, technology, and engineering, ranging from gravity to mechanics. In general, non-linear PDEs are important tools that can be used in various fields such as plasma physics, mathematical biology, solid state physics, and fluid dynamics for modeling nonlinear dynamic phenomena [5]. The majority of dynamic schemes can be denoted by an acceptable array of PDEs. It is also well-appreciated that PDEs, such as Poincare and Calabi conjecture models, are utilized to solve mathematical difficulties.

It has been found that the non-linear development of shallow water waves in the fluid dynamics is described by utilizing the coupled scheme Whitham–Broer–Kaup equations (WBKEs) [6]. The coupled scheme of the above equations was developed by Whitham, Broer, and Kaup [7–9]. The above equation defines the propagation of shallow water waves with specific diffusion families.

In the classical order, the major equations of the said phenomena are given as:

$$\begin{aligned} D^\delta\_\eta \mu(a,\eta) + \mu(a,\eta) \frac{\partial \mu(a,\eta)}{\partial a} + \frac{\partial \mu(a,\eta)}{\partial a} + g \frac{\partial \nu(a,\eta)}{\partial a} &= 0, \\ D^\delta\_\eta \nu(a,\eta) + \mu(a,\eta) \frac{\partial \nu(a,\eta)}{\partial a} + \nu(a,\eta) \frac{\partial \mu(a,\eta)}{\partial a} + p \frac{\partial^3 \mu(a,\eta)}{\partial a^3} - g \frac{\partial^2 \nu(a,\eta)}{\partial a^2} &= 0, \quad 0 < \delta \le 1. \end{aligned}$$

Here *μ*(*α*, *η*) and *ν*(*α*, *η*) describe the straight velocity and height, which deviate from the equilibrium situation of the fluid, respectively, and *p* and *q* are constants expressed in various diffusion forces. Investigating solutions to such nonlinear PDEs over the last several decades it is an important research area [10]. Several scientists have developed numerous mathematical techniques to explore the approximate solutions to nonlinear PDEs. Aminikhah and Biazar [11] used the HPM (homotopy perturbation method) to solve the coupled model of Brusselator and Burger equations. Noor and Mohyud-Din [12] utilized HPM to examine the solutions of different classical orders of PDEs. Ahmad et al. [5] studied a coupled scheme result of WBKEs by the Adomian decomposition method (ADM). Whitham–Broer–Kaup equations are solved by other researchers using different analytical and numerical methods, such as the hyperbolic function method [13], residual power series method (RPSM) [14], Adomian decomposition method [15], reduced differential transformation method [16], homotopy perturbation method [17,18], exp-function method [19], Lie Symmetry analysis [20,21], G /G2-Expansion method [22], and homotopy analysis method [23]. Recently, Amjad et al. [10] used the result of a standard order coupled of fractional-order Whitham–Broer–Kaup equation by the Laplace decomposition method.

Singh et al. [24] suggested the q-HATM, which is a well-designed mixture of Laplace transform and q-HAM. The future system monitors and manipulates the sequences result, which converges quickly to the exact solutions for the problem. The strength of the proposed technique is its ability to combine two powerful algorithms to solve both numerically and analytically linear and non-linear fractional-order differential equations. A future procedure has several study properties that include a non-local effect, straightforward result system, promising broad convergence area, and free of any perturbation, discretization and assumption. It is worth disclosing that, using semi-analytical techniques, the Laplace transform takes less C.P.U. time to determine the solutions of complex nonlinear models and phenomena that occur in technology and science. The solution q-HATM includes two auxiliary parameters *h* and *n*, which aims to help us modify and control the solution's convergence [25]. Recently, with the help of q-HATM, several researchers studied different phenomena in different fields for example, Singh et al. studied to find the advection-dispersion equation solution [26] and Srivastava et al. used an arbitrary order vibration equation model [27].

The natural decomposition method (NDM) is a mixture of the Adomian decomposition method and the natural transform method (NTM). In 2014, S. Maitama and M. Rawashdeh first implemented the NDM [28,29] to solve linear and non-linear ordinary differential equations (ODEs) and PDEs that occur in several fields of science. A huge quantity of physical models have been studied using NDM, such as the study of fractional order diffusion equations [30], fractional order delay PDEs [31] nonlinear PDEs [32,33], the fractional uncertain flow of a system of polytropic gas [34], fractional-order physical schemes [35], fractional wave and heat problems [36], and fractional telegraph equation [37].

In the current research article, two analytical methods, namely the natural decomposition method and q-homotopy analysis transform method are used to solve the fractional-order Whitham–Broer–Kaup equation. The solutions obtained by the proposed techniques are very simple and straightforward. Moreover, the accuracy of the present methods is sufficient to obtain the analytical solution of the targeted problems. The obtained solutions are compared and to found to be in a good agreement with the exact solution for the problem. This article introduces an approximate analytical solution of a multi dimensional, time fractional model of the Whitham–Broer–Kaup equation by implementing NDM and q-HATM.

#### **2. Preliminaries Concepts**

**Definition 1.** *The Laplace transformation of a Caputo fractional derivative Dδg*(*η*) *is described as:*

$$\mathcal{E}[D^{\delta}g(\eta)] = s^{\delta}\mathcal{R}(s) - \sum\_{j=0}^{m-1} s^{\delta - (j+1)} [g^{(j)}(0^{+})] \quad m-1 \le \delta < m$$

**Definition 2.** *The natural transformation of the <sup>g</sup>*(*η*) *function is represented by <sup>N</sup>*+[*g*(*η*)] *for <sup>η</sup>* <sup>∈</sup> *<sup>R</sup> and identified by:*

$$N^{+}[\mathcal{g}(\eta)] = \mathcal{R}(\mathfrak{s}, \mathfrak{u}) = \int\_{-\infty}^{\infty} e^{-s\eta} \mathcal{g}(\eta) d\eta; \quad \mathfrak{s}, \mathfrak{u} \in (-\infty, \infty),$$

*where the natural transformation variables are s and u. If g*(*η*)*Q*(*η*) *is described on the real positive axis, the natural transformation is described as:*

$$N^+[\mathcal{g}(\eta)Q(\eta)] = N^+[\mathcal{g}(\eta)] = \mathcal{R}^+(s, u) = \int\_0^\infty e^{-s\eta} \mathcal{g}(\eta) d\eta; \quad s, u \in (0, \infty), \quad and \quad \eta \in \mathcal{R}$$

*where Q*(*η*) *represents the function of Heaviside. Simply, for u* = 1*, the equation is reduced to the Laplace transformation, and for s* = 1*, the equation is the Sumud transformation.*

**Theorem 1.** *Let* R(*s*, *u*) *be the natural transformation of the function g*(*η*)*, then the natural transform* <sup>R</sup>*δ*(*s*, *<sup>u</sup>*) *of the Riemann–Liouville fractional derivative of g*(*η*) *is symbolized by Dδg*(*η*) *and is presented as:*

$$N^{+}[D^{\delta}g(\eta)] = \mathcal{R}\_{\delta}(\mathbf{s},\boldsymbol{\mu}) = \frac{\mathbf{s}^{\delta}}{\mathbf{u}^{\delta}}\mathcal{R}(\mathbf{s},\boldsymbol{\mu}) - \sum\_{j=0}^{m-1} \frac{\mathbf{s}^{j}}{\mathbf{u}^{\delta-j}}[D^{\delta-j-1}g(\eta)]\_{\eta=0,\boldsymbol{\mu}}$$

*where δ is the order and m be any positive integer. Furthermore, m* − 1 ≤ *δ* < *m.*

**Theorem 2.** *Let* R(*s*, *u*) *be the natural transformation of the g*(*η*)*, then the natural transformation* R*δ*(*s*, *u*) *of the Caputo fractional derivative of g*(*η*) *is symbolized by Dδg*(*η*) *and is represented as:*

$$N^{+}[^{c}D^{\delta}g(\eta)] = \mathcal{R}^{c}\_{\delta}(s,\mu) = \frac{s^{\delta}}{\iota^{\delta}}\mathcal{R}(s,\mu) - \sum\_{j=0}^{m-1} \frac{s^{\delta-(j+1)}}{\iota^{\delta-j}}[D^{j}g(\eta)]\_{\eta=0} \quad m-1 \le \delta < m$$

**Definition 3.** *The fractional derivative of g* <sup>∈</sup> *<sup>C</sup><sup>m</sup>* <sup>−</sup><sup>1</sup> *in the Caputo sense is represented as:*

$$D\_{\eta}^{\delta}g(\eta) = \begin{cases} \frac{\partial^{n}\varrho(\eta)}{\partial\eta^{n}}, & \delta = m \in N\_{\prime} \\ \frac{1}{\Gamma(m-\delta)} \int\_{0}^{\eta} (\eta-\phi)^{m-\delta-1} g^{m}(\phi) \partial\phi, & m-1 < \delta < m, \quad m \in N. \end{cases}$$

**Definition 4.** *Function of Mittag–Leffler, Eδ*(*b*) *for δ* > 0 *is defined as:*

$$E\_{\delta}(b) = \sum\_{m=0}^{\infty} \frac{b^m}{\Gamma(\delta m + 1)} \qquad \delta > 0 \quad b \in \mathbb{C}\_{\prime}$$

#### **3. The Procedure of NDM**

In this section, we describe the NDM solution scheme for fractional partial differential equations.

$$\begin{aligned} D\_{\eta}^{\delta}\mu(\mathfrak{a},\eta) + \mathcal{R}\_{1}(\mu,\nu) + \mathcal{N}\_{1}(\mu,\nu) - \mathcal{P}\_{1}(\mathfrak{a},\eta) &= 0, \\ D\_{\eta}^{\delta}\nu(\mathfrak{a},\eta) + \mathcal{R}\_{2}(\mu,\nu) + \mathcal{N}\_{2}(\mu,\nu) - \mathcal{P}\_{2}(\mathfrak{a},\eta) &= 0, \quad 0 < \delta \le 1, \end{aligned} \tag{1}$$

with the initial condition:

$$
\mu(a,0) = \mathfrak{g}\_1(a), \quad \nu(a,0) = \mathfrak{g}\_2(a). \tag{2}
$$

where is *D<sup>δ</sup> <sup>η</sup>* <sup>=</sup> *<sup>∂</sup><sup>δ</sup> ∂η<sup>δ</sup>* the Caputo fractional derivative of order *δ*, R1, R<sup>2</sup> and N1, N<sup>2</sup> are linear and non-linear functions or operators, respectively, and P1,P<sup>2</sup> are source functions.

Applying the natural transform to Equation (1),

$$\begin{split} &N^{+}\left[D^{\delta}\_{\eta}\mu(\mathfrak{a},\eta)\right] + N^{+}\left[\mathcal{R}\_{1}(\mu,\nu) + \mathcal{N}\_{1}(\mu,\nu) - \mathcal{P}\_{1}(\mathfrak{a},\eta)\right] = 0, \\ &N^{+}\left[D^{\delta}\_{\eta}\nu(\mathfrak{a},\eta)\right] + N^{+}\left[\mathcal{R}\_{2}(\mu,\nu) + \mathcal{N}\_{2}(\mu,\nu) - \mathcal{P}\_{2}(\mathfrak{a},\eta)\right] = 0. \end{split} \tag{3}$$

Using the differentiation property of natural transform, we get:

$$\begin{split} \left[N^{+}\left[\mu(\mathfrak{a},\eta)\right]\right] &= \frac{u^{\delta}}{\mathfrak{s}^{\delta}}\sum\_{k=0}^{m-1} \frac{\mathfrak{s}^{\delta-k-1}}{\mathfrak{u}^{\delta-k}} \frac{\partial^{k}\mu(\mathfrak{a},\eta)}{\partial^{k}\eta}|\_{\eta=0} + \frac{u^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left[\mathcal{P}\_{1}(\mathfrak{a},\eta)\right] - \frac{u^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left\{\mathcal{R}\_{1}(\mu,\nu) + \mathcal{N}\_{1}(\mu,\nu)\right\}|\_{\nu} \\ \left[N^{+}\left[\nu(\mathfrak{a},\eta)\right]\right] &= \frac{u^{\delta}}{\mathfrak{s}^{\delta}}\sum\_{k=0}^{m-1} \frac{\mathfrak{s}^{\delta-k-1}}{\mathfrak{u}^{\delta-k}} \frac{\partial^{k}\nu(\mathfrak{a},\eta)}{\partial^{k}\eta}|\_{\eta=0} + \frac{u^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left[\mathcal{P}\_{2}(\mathfrak{a},\eta)\right] - \frac{u^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left\{\mathcal{R}\_{2}(\mu,\nu) + \mathcal{N}\_{2}(\mu,\nu)\right\}|\_{\nu} \end{split} (4)$$

NDM describes the solution of infinite series *μ*(*α*, *η*) and *ν*(*α*, *η*),

$$\mu(a,\eta) = \sum\_{m=0}^{\infty} \mu\_{\text{ll}}(a,\eta), \quad \nu(a,\eta) = \sum\_{m=0}^{\infty} \nu\_{\text{ll}}(a,\eta), \tag{5}$$

Adomian polynomials of non-linear terms of N<sup>1</sup> and N<sup>2</sup> are represented as:

$$\mathcal{N}\_1(\mu, \nu) = \sum\_{m=0}^{\infty} \mathcal{A}\_{m\nu} \quad \mathcal{N}\_2(\mu, \nu) = \sum\_{m=0}^{\infty} \mathcal{B}\_{m\nu} \tag{6}$$

All forms of non-linearity of the Adomian polynomials can be defined as:

$$\begin{split} \mathcal{A}\_{m} &= \frac{1}{m!} \left[ \frac{\partial^{m}}{\partial \lambda^{m}} \left\{ \mathcal{N}\_{1} \left( \sum\_{k=0}^{\infty} \lambda^{k} \mu\_{k\prime} \sum\_{k=0}^{\infty} \lambda^{k} \nu\_{k} \right) \right\} \right]\_{\lambda=0} \\ \mathcal{B}\_{m} &= \frac{1}{m!} \left[ \frac{\partial^{m}}{\partial \lambda^{m}} \left\{ \mathcal{N}\_{2} \left( \sum\_{k=0}^{\infty} \lambda^{k} \mu\_{k\prime} \sum\_{k=0}^{\infty} \lambda^{k} \nu\_{k} \right) \right\} \right]\_{\lambda=0} \end{split} \tag{7}$$

Substituting Equations (13) and (14) into Equation (12) gives:

$$\begin{split} N^{+}\left[\sum\_{m=0}^{\infty}\mu\_{\rm{tr}}(a,\eta)\right] &= \frac{\mu^{\delta}}{\mathfrak{s}^{\delta}}\sum\_{k=0}^{m-1}\frac{\mathfrak{s}^{\delta-k-1}}{\mathfrak{u}^{\delta-k}}\frac{\partial^{\delta}\mu(a,\eta)}{\partial^{\delta}\eta}|\_{\eta=0} + \frac{\mathfrak{u}^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left\{\mathcal{P}\_{\rm{T}}(a,\eta)\right\} - \frac{\mathfrak{u}^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left\{\mathcal{R}\_{\rm{L}}\left(\sum\_{m=0}^{\infty}\mu\_{\rm{tr}}\sum\_{m=0}^{\infty}\nu\_{m}\right) + \sum\_{m=0}^{\infty}\mathcal{A}\_{\rm{m}}\right\}, \\ N^{+}\left[\sum\_{m=0}^{\infty}\nu\_{\rm{tr}}(a,\eta)\right] &= \frac{\mathfrak{u}^{\delta}}{\mathfrak{s}^{\delta}}\sum\_{k=0}^{m-1}\frac{\mathfrak{s}^{\delta-k-1}}{\mathfrak{u}^{\delta-k}}\frac{\partial^{\delta}\nu(a,\eta)}{\partial^{\delta}\eta}|\_{\eta=0} + \frac{\mathfrak{u}^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left\{\mathcal{P}\_{\rm{T}}(a,\eta)\right\} - \frac{\mathfrak{u}^{\delta}}{\mathfrak{s}^{\delta}}N^{+}\left\{\mathcal{R}\_{\rm{L}}\left(\sum\_{m=0}^{\infty}\mu\_{m},\sum\_{m=0}^{\infty}\nu\_{m}\right) + \sum\_{m=0}^{\infty}\mathcal{B}\_{\rm{m}}\right\}, \end{split} \tag{8}$$

Applying the inverse natural transformation of Equation (16),

$$\begin{split} \sum\_{m=0}^{\infty} \mu\_{m}(\boldsymbol{a}, \boldsymbol{\eta}) &= N^{-} \left[ \frac{\mu^{\delta}}{s^{\delta}} \sum\_{k=0}^{m-1} \frac{\delta^{\delta-k-1}}{\mu^{\delta-k}} \frac{\partial^{k} \mu(\boldsymbol{a}, \boldsymbol{\eta})}{\partial^{k} \boldsymbol{\eta}} \vert\_{\boldsymbol{\eta}=0} + \frac{\mu^{\delta}}{s^{\delta}} N^{+} \{\mathcal{P}\_{1}(\boldsymbol{a}, \boldsymbol{\eta})\} \right] - N^{-} \left[ \frac{\mu^{\delta}}{s^{\delta}} N^{+} \{\mathcal{R}\_{1} \left( \sum\_{m=0}^{\infty} \mu\_{m} \sum\_{m=0}^{\infty} \boldsymbol{\nu}\_{m} \right) + \sum\_{m=0}^{\infty} \mathcal{A}\_{m} \right] \,, \\ \sum\_{m=0}^{\infty} \boldsymbol{\nu}\_{m}(\boldsymbol{a}, \boldsymbol{\eta}) &= N^{-} \left[ \frac{\mu^{\delta}}{s^{\delta}} \sum\_{k=0}^{m-1} \frac{\mathcal{S}^{\delta-k-1}}{\mu^{\delta-k}} \frac{\partial^{k} \nu(\boldsymbol{a}, \boldsymbol{\eta})}{\partial^{k} \boldsymbol{\eta}} \vert\_{\boldsymbol{\eta}=0} + \frac{\mu^{\delta}}{s^{\delta}} N^{+} \{\mathcal{P}\_{2}(\boldsymbol{a}, \boldsymbol{\eta})\} \right] - N^{-} \left[ \frac{\mu^{\delta}}{s^{\delta}} N^{+} \{\mathcal{R}\_{2} \left( \sum\_{m=0}^{\infty} \mu\_{m} \sum\_{m=0}^{\infty} \boldsymbol{\nu}\_{m} \right) + \sum\_{m=0}^{\infty} \mathcal{B}\_{m} \right] \,, \end{split} \tag{9}$$

we define the following terms,

*μ*0(*α*, *η*) = *N*−[ *uδ sδ m*−1 ∑ *k*=0 *sδ*−*k*−<sup>1</sup> *uδ*−*<sup>k</sup> ∂kμ*(*α*, *η*) *<sup>∂</sup>k<sup>η</sup>* <sup>|</sup>*η*=<sup>0</sup> <sup>+</sup> *<sup>u</sup><sup>δ</sup> <sup>s</sup><sup>δ</sup> <sup>N</sup>*+{P1(*α*, *<sup>η</sup>*)}], *ν*0(*α*, *η*) = *N*−[ *uδ sδ m*−1 ∑ *k*=0 *sδ*−*k*−<sup>1</sup> *uδ*−*<sup>k</sup> ∂kν*(*α*, *η*) *<sup>∂</sup>k<sup>η</sup>* <sup>|</sup>*η*=<sup>0</sup> <sup>+</sup> *<sup>u</sup><sup>δ</sup> <sup>s</sup><sup>δ</sup> <sup>N</sup>*+{P2(*α*, *<sup>η</sup>*)}], (10) *μ*1(*α*, *η*) = −*N*−[ *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*+{R1(*μ*0, *<sup>ν</sup>*0) + <sup>A</sup>0}], *ν*1(*α*, *η*) = −*N*−[ *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*+{R2(*μ*0, *<sup>ν</sup>*0) + <sup>B</sup>0}],

the general for *m* ≥ 1, is given by:

$$\begin{aligned} \mu\_{m+1}(\mathfrak{a},\eta) &= -N^- [\frac{\mathfrak{u}^\delta}{\mathfrak{s}^\delta} N^+ \{ \mathcal{R}\_1(\mu\_{m'} \nu\_{m}) + \mathcal{A}\_m \} ], \\ \nu\_{m+1}(\mathfrak{a},\eta) &= -N^- [\frac{\mathfrak{u}^\delta}{\mathfrak{s}^\delta} N^+ \{ \mathcal{R}\_2(\mu\_{m'} \nu\_{m}) + \mathcal{B}\_m \} ]. \end{aligned}$$

#### **4. Fundamental Idea of q-Homotopy Analysis Transform Method**

To introduce the basic concept of the current method, we consider a fractional-order nonlinear PDEs of the form:

$$D\_{\eta}^{\delta}\mu(\mathfrak{a},\mathfrak{E},\eta) + R\mu(\mathfrak{a},\mathfrak{E},\eta) + N\mu(\mathfrak{a},\mathfrak{E},\eta) = f(\mathfrak{a},\mathfrak{E},\eta), \quad 1 < \delta \le n,\tag{11}$$

where *D<sup>δ</sup> ημ*(*α*, *β*, *η*) denote's the Caputo's fractional derivative *R* and *N* are linear and non-linear functions or operators. Using the differentiation property of the Laplace transform on Equation (12), we get:

$$\mathrm{s}^{\delta}\mathcal{L}[\mu(\mathfrak{a},\mathfrak{\beta},\mathfrak{\eta})] - \sum\_{k=0}^{m-1} \mathrm{s}^{\delta-k-1} \frac{\mathrm{\hat{\sigma}}^{k}\mu(\mathfrak{a},\mathfrak{\beta},\mathfrak{\eta})}{\mathrm{\hat{\sigma}}^{k}\mathfrak{\eta}}|\_{\mathfrak{\eta}=0} + \mathfrak{e}[\mathrm{R}\mu(\mathfrak{a},\mathfrak{\beta},\mathfrak{\eta}) + \mathrm{N}\mu(\mathfrak{a},\mathfrak{\beta},\mathfrak{\eta})] = \mathfrak{e}[f(\mathfrak{a},\mathfrak{\beta},\mathfrak{\eta})],\tag{12}$$

$$\mathbb{E}[D\_{\eta}^{\delta}\mu(\mathfrak{a},\mathfrak{b},\eta)] = \mathbb{E}[\mu(\mathfrak{a},\mathfrak{b},\eta)] - \frac{1}{\mathfrak{s}^{\delta}} \sum\_{k=0}^{m-1} \mathfrak{s}^{\delta-k-1} \frac{\partial^{k}\mu(\mathfrak{a},\mathfrak{b},\eta)}{\partial^{k}\eta}|\_{\eta=0}.\tag{13}$$

On simplifying Equation (13), we have:

$$\left[\mathcal{L}\mu(\boldsymbol{a},\boldsymbol{\beta},\boldsymbol{\eta})-\frac{1}{s^{\delta}}\sum\_{k=0}^{m-1}s^{\delta-k-1}\frac{\partial^{k}\mu(\boldsymbol{a},\boldsymbol{\beta},\boldsymbol{\eta})}{\partial^{k}\boldsymbol{\eta}}\vert\_{\boldsymbol{\eta}=\boldsymbol{0}}+\frac{1}{s^{\delta}}\mathcal{E}[\boldsymbol{R}\mu(\boldsymbol{a},\boldsymbol{\beta},\boldsymbol{\eta})+\mathcal{N}\mu(\boldsymbol{a},\boldsymbol{\beta},\boldsymbol{\eta})-f(\boldsymbol{a},\boldsymbol{\beta},\boldsymbol{\eta})] \right] = 0. \tag{14}$$

We can describe the non-linear operator as:

$$\begin{split} N[\phi(a,\boldsymbol{\beta},\boldsymbol{\eta};q)] &= \mathfrak{E}[\phi(a,\boldsymbol{\beta},\boldsymbol{\eta};q)] - \frac{1}{s^{\delta}} \sum\_{k=0}^{m-1} s^{\delta-k-1} \frac{\partial^{k} \phi(a,\boldsymbol{\beta},\boldsymbol{\eta};q)}{\partial^{k} \boldsymbol{\eta}}|\_{\boldsymbol{\eta}=0} + \frac{1}{s^{\delta}} \mathfrak{E}[\mathcal{R}\phi(a,\boldsymbol{\beta},\boldsymbol{\eta};q)] \\ &+ \frac{1}{s^{\delta}} \mathfrak{E}[\mathcal{N}\phi(a,\boldsymbol{\beta},\boldsymbol{\eta};q)] - \frac{1}{s^{\delta}} \mathfrak{E}[f(a,\boldsymbol{\beta},\boldsymbol{\eta})], \end{split} \tag{15}$$

where *<sup>q</sup>* <sup>∈</sup> [0, <sup>1</sup> *<sup>n</sup>* ], and *φ*(*α*, *β*, *η*; *q*) is real function of *α*, *β*, *η*, and *q*. The concept of a nonzero auxiliary function of homotopy is the following:

$$(1 - nq)\mathcal{E}[\phi(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}; q) - \mu\_0(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) = \hbar q H(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) N[\phi(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}, q)],\tag{16}$$

where *£* a sign of the Laplace transformation, *<sup>q</sup>* <sup>∈</sup> [0, <sup>1</sup> *<sup>n</sup>* ](*n* ≥ 1) is the embedding parameter, *h*¯ = 0 is an auxiliary parameter, *H*(*α*, *β*, *η*) signifies a nonzero auxiliary function, *φ*(*α*, *β*, *η*; *q*) is an unidentified function, and *μ*0(*α*, *β*, *η*) is an initial guess of *μ*(*α*, *β*, *η*). The subsequent outcomes hold correspondingly for *q* = 0 and *q* = <sup>1</sup> *n* .

$$
\phi(a,\beta,\eta;0) = \mu\_0(a,\beta,\eta), \\
\phi(a,\beta,\eta,\frac{1}{n}) = \mu(a,\beta,\eta). \tag{17}
$$

Thus, by intensifying *q* from 0 to <sup>1</sup> *<sup>n</sup>* , the result *φ*(*α*, *β*, *η*; *q*) converge from *μ*0(*α*, *β*, *η*) to the solution *μ*(*α*, *β*, *η*). Expand the function *φ*(*α*, *β*, *η*, *q*) in sequences form by using the Taylor theorem near to *q*, where one can get:

$$\phi(a,\beta,\eta;\eta) = \mu\_0(a,\beta,\eta) + \sum\_{m=1}^{\infty} \mu\_m(a,\beta,\eta)q^m,\tag{18}$$

where,

$$\mu\_0(\mathfrak{a}, \mathfrak{E}, \eta) = \frac{1}{m!} \frac{\partial^m \phi(\mathfrak{a}, \mathfrak{E}, \eta; q)}{\partial q^m} \mid q = 0 \tag{19}$$

On selecting the auxiliary linear operator, *μ*0(*α*, *β*, *η*) , n and *h*¯, the series (19) converge at *q* = <sup>1</sup> *<sup>n</sup>* and then it produces one of the results for Equation (12):

$$
\mu(a,\beta,\eta) = \mu\_0(a,\beta,\eta) + \sum\_{m=1}^{\delta} \mu\_m(a,\beta,\eta) (\frac{1}{n})^m \tag{20}
$$

Now, differentiating the zero-th order distortion Equation (17) m-times with respect to *q* and then dividing by *m*! and lastly taking *q* = 0, which provides:

$$\mathcal{E}[\mu\_m(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) - \mathcal{K}\_m \mu\_{m-1}(\mathfrak{a}, \mathfrak{z}, \mathfrak{y})] = \hbar \mathfrak{R}\_m(\mu\_{m-1}^\rightarrow),\tag{21}$$

where,

$$
\mu\_m^{\rightarrow} = \mu\_0(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) + \mu\_1(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) \,. \tag{22}
\\
\mu\_m(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) . \tag{22}
$$

Using the inverse Laplace transformation on Equation (22), it produces:

$$
\mu\_m(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) = \mathsf{K}\_m \mu\_{m-1}(\mathfrak{a}, \mathfrak{z}, \mathfrak{y}) + \hbar \mathsf{L}^{-1} [\mathfrak{R}\_m(\mu\_{m-1}^\rightarrow)] \tag{23}
$$

where,

$$\begin{split} \mathbb{R}(\mu\_{m-1}^{\rightarrow}) &= \mathbb{\varepsilon}[\mu\_{m-1}(a,\beta,\eta)] - (1 - \frac{K\_{\text{m}}}{n})(\sum\_{k=0}^{m-1} s^{\delta-k-1} \frac{\partial^k \mu(a,\beta,\eta)}{\partial^k \eta}|\_{\eta=0} + \frac{1}{s^\beta} \mathbb{\varepsilon}[f(a,\beta,\eta)]) \\ &+ \frac{1}{s^\beta} \mathbb{\varepsilon}[\mathbb{R}(\mu\_{m-1} + H\_{m-1}), \end{split} \tag{24}$$

And,

$$k\_m = \begin{cases} 0, m \le 1 \\ 1, m > 1 \end{cases} \tag{25}$$

In Equation (25), *Hm* denotes a homotopy polynomial and is defined as:

$$H\_m = \mu\_0(a, \beta, \eta) = \frac{1}{m!} \frac{\partial^m \phi(a, \beta, \eta; q)}{\partial q^m} \mid q = 0 \quad \text{and} \quad \phi(a, \beta, \eta; q) = \phi\_0 + q\phi\_1 + q^2\phi\_2 + \dots \dots \tag{26}$$

By Equations (24) and (25), we have:

$$\begin{split} \mu\_{m}(a,\boldsymbol{\beta},\boldsymbol{\eta}) &= (k\_{m} + \hbar)\mu\_{m-1}(a,\boldsymbol{\beta},\boldsymbol{\eta}) - (1 - \frac{k\_{m}}{n})\mathcal{L}^{-1}[\sum\_{k=0}^{m-1} s^{\delta-k-1} \frac{\partial^{k}\mu(a,\boldsymbol{\beta},\boldsymbol{\eta})}{\partial^{k}\boldsymbol{\eta}}|\_{\boldsymbol{\eta}=0} + \frac{1}{s^{\delta}}\mathcal{L}[f(a,\boldsymbol{\beta},\boldsymbol{\eta})]) \\ &+ \hbar\mathcal{L}^{-1}[\frac{1}{s^{\delta}}\mathcal{L}[\mathcal{R}(\mu\_{m-1} + H\_{m-1})]. \end{split} \tag{27}$$

On solving Equation (28) for *m* = 1, 2, 3, 4, ...... with the help of *μ*0(*α*, *β*, *η*) = *μ*(*x*, *y*, 0) and Equation (25), we get the iterative terms of *μm*(*α*, *β*, *η*). The q-homotopy analysis transform method series solution is given by:

$$
\mu(\alpha, \beta, \eta) = \sum\_{m=0}^{\infty} \mu\_m(\alpha, \beta, \eta). \tag{28}
$$

#### **5. Numerical Examples**

**Example 1.** *Consider the coupled system of the fractional-order Whitham–Broer–Kaup equations with:*

$$\begin{cases} D^{\delta}\_{\eta}\mu(a,\eta) + \mu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} + \frac{\partial\mu(a,\eta)}{\partial a} + \frac{\partial\nu(a,\eta)}{\partial a} = 0, \\ D^{\delta}\_{\eta}\nu(a,\eta) + \mu(a,\eta)\frac{\partial\nu(a,\eta)}{\partial a} + \nu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} + 3\frac{\partial^{3}\mu(a,\eta)}{\partial a^{3}} - \frac{\partial^{2}\nu(a,\eta)}{\partial a^{2}} = 0, \\ 0 < \delta \le 1, \quad -1 < \eta \le 1, \quad -10 \le a \le 10, \end{cases} \tag{29}$$

*with the initial condition:*

$$\begin{cases} \quad \mu(\alpha, 0) = \frac{1}{2} - 8 \tanh(-2\alpha), \\\quad \nu(\alpha, 0) = 16 - 16 \tanh^2(-2\alpha). \end{cases} \tag{30}$$

*Firstly , we will solve this scheme by using the NDM.*

After the natural transformation of Equation (29), we get:

*N*<sup>+</sup> *∂δμ*(*α*, *η*) *∂η<sup>δ</sup>* <sup>=</sup> <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂μ*(*α*, *η*) *∂α* <sup>+</sup> *∂μ*(*α*, *<sup>η</sup>*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* " , *N*<sup>+</sup> *∂δν*(*α*, *η*) *∂η<sup>δ</sup>* <sup>=</sup> <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂ν*(*α*, *η*) *∂α* <sup>+</sup> *<sup>ν</sup>*(*α*, *<sup>η</sup>*) *∂μ*(*α*, *η*) *∂α* <sup>+</sup> <sup>3</sup> *∂*3*μ*(*α*, *η*) *∂α*<sup>3</sup> <sup>−</sup> *<sup>∂</sup>*2*ν*(*α*, *<sup>η</sup>*) *∂α*<sup>2</sup> " , *sδ <sup>u</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> {*μ*(*α*, *<sup>η</sup>*)} <sup>−</sup> *<sup>s</sup>δ*−<sup>1</sup> *<sup>u</sup><sup>δ</sup> <sup>μ</sup>*(*α*, 0) = <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂μ*(*α*, *η*) *∂α* <sup>+</sup> *∂μ*(*α*, *<sup>η</sup>*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* " *sδ <sup>u</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> {*ν*(*α*, *<sup>η</sup>*)} <sup>−</sup> *<sup>s</sup>δ*−<sup>1</sup> *<sup>u</sup><sup>δ</sup> <sup>ν</sup>*(*α*, 0) = <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂ν*(*α*, *η*) *∂α* <sup>+</sup> *<sup>ν</sup>*(*α*, *<sup>η</sup>*) *∂μ*(*α*, *η*) *∂α* <sup>+</sup> <sup>3</sup> *∂*3*μ*(*α*, *η*) *∂α*<sup>3</sup> <sup>−</sup> *<sup>∂</sup>*2*ν*(*α*, *<sup>η</sup>*) *∂α*<sup>2</sup> " ,

The above algorithm is reduced to be simplified:

$$\begin{split} \{N^{+}\left\{\mu(a,\eta)\right\}\} &= \frac{1}{s} \left\{\mu(a,0)\right\} - \frac{\mathrm{n}^{\delta}}{\mathrm{s}^{\delta}} N^{+} \left[\mu(a,\eta) \frac{\partial \mu(a,\eta)}{\partial a} + \frac{\partial \mu(a,\eta)}{\partial a} + \frac{\partial \nu(a,\eta)}{\partial a}\right], \\ \{N^{+}\left\{\nu(a,\eta)\right\}\} &= \frac{1}{s} \left\{\nu(a,0)\right\} - \frac{\mathrm{n}^{\delta}}{\mathrm{s}^{\delta}} N^{+} \left[\mu(a,\eta) \frac{\partial \nu(a,\eta)}{\partial a} + \nu(a,\eta) \frac{\partial \mu(a,\eta)}{\partial a} + 3 \frac{\partial^{3} \mu(a,\eta)}{\partial a^{3}} - \frac{\partial^{2} \nu(a,\eta)}{\partial a^{2}}\right], \end{split} \tag{31}$$

Applying inverse natural transformation, we get:

$$\begin{split} \mu(\boldsymbol{a},\boldsymbol{\eta}) &= \mu(\boldsymbol{a},\boldsymbol{0}) - \boldsymbol{N}^{-} \left[ \frac{\boldsymbol{\mu}^{\delta}}{s^{\delta}} \boldsymbol{N}^{+} \left[ \mu(\boldsymbol{a},\boldsymbol{\eta}) \frac{\partial \mu(\boldsymbol{a},\boldsymbol{\eta})}{\partial \boldsymbol{a}} + \frac{\partial \mu(\boldsymbol{a},\boldsymbol{\eta})}{\partial \boldsymbol{a}} + \frac{\partial \nu(\boldsymbol{a},\boldsymbol{\eta})}{\partial \boldsymbol{a}} \right] \right], \\ \nu(\boldsymbol{a},\boldsymbol{\eta}) &= \nu(\boldsymbol{a},\boldsymbol{0}) - \boldsymbol{N}^{-} \left[ \frac{\boldsymbol{\mu}^{\delta}}{s^{\delta}} \boldsymbol{N}^{+} \left[ \mu(\boldsymbol{a},\boldsymbol{\eta}) \frac{\partial \nu(\boldsymbol{a},\boldsymbol{\eta})}{\partial \boldsymbol{a}} + \nu(\boldsymbol{a},\boldsymbol{\eta}) \frac{\partial \mu(\boldsymbol{a},\boldsymbol{\eta})}{\partial \boldsymbol{a}} + 3 \frac{\partial^{3} \mu(\boldsymbol{a},\boldsymbol{\eta})}{\partial \boldsymbol{a}^{3}} - \frac{\partial^{2} \nu(\boldsymbol{a},\boldsymbol{\eta})}{\partial \boldsymbol{a}^{2}} \right] \right], \end{split} \tag{32}$$

Assume that the unknown functions *μ*(*α*, *η*) and *ν*(*α*, *η*) infinite series solution is as follows:

$$\mu(\mathfrak{a},\eta) = \sum\_{m=0}^{\infty} \mu\_m(\mathfrak{a},\eta), \text{ and } \quad \nu(\mathfrak{a},\eta) = \sum\_{m=0}^{\infty} \nu\_m(\mathfrak{a},\eta).$$

Remember that *μμα* = ∑<sup>∞</sup> *<sup>m</sup>*=<sup>0</sup> <sup>A</sup>*m*, *μνα* <sup>=</sup> <sup>∑</sup><sup>∞</sup> *<sup>m</sup>*=<sup>0</sup> <sup>B</sup>*<sup>m</sup>* and *νμα* <sup>=</sup> <sup>∑</sup><sup>∞</sup> *<sup>m</sup>*=<sup>0</sup> C*<sup>m</sup>* are the Adomian polynomials and the nonlinear terms were characterized. Using such terms, Equation (32) can be rewritten in the form:

∞ ∑ *m*=0 *μm*(*α*, *η*) = *μ*(*α*, 0) − *N*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 <sup>A</sup>*<sup>m</sup>* <sup>+</sup> *∂μ*(*α*, *<sup>η</sup>*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* , ∞ ∑ *m*=0 *νm*(*α*, *η*) = *ν*(*α*, 0) − *N*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 B*<sup>m</sup>* + ∞ ∑ *m*=0 C*<sup>m</sup>* + 3 *∂*3*μ*(*α*, *η*) *∂α*<sup>3</sup> <sup>−</sup> *<sup>∂</sup>*2*ν*(*α*, *<sup>η</sup>*) *∂α*<sup>2</sup> , ∞ ∑ *m*=0 *<sup>μ</sup>m*(*α*, *<sup>η</sup>*) = <sup>1</sup> <sup>2</sup> <sup>−</sup> 8 tanh(−2*α*) <sup>−</sup> *<sup>N</sup>*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 <sup>A</sup>*<sup>m</sup>* <sup>+</sup> *∂μ*(*α*, *<sup>η</sup>*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* , ∞ ∑ *m*=0 *<sup>ν</sup>m*(*α*, *<sup>η</sup>*) = <sup>16</sup> <sup>−</sup> 16 tanh2(−2*α*) <sup>−</sup> *<sup>N</sup>*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 B*<sup>m</sup>* + ∞ ∑ *m*=0 C*<sup>m</sup>* + 3 *∂*3*μ*(*α*, *η*) *∂α*<sup>3</sup> <sup>−</sup> *<sup>∂</sup>*2*ν*(*α*, *<sup>η</sup>*) *∂α*<sup>2</sup> , (33)

According to Equation (7), all forms of non-linearity the Adomian polynomials can be defined as:

$$\begin{split} \mathcal{A}\_{0} &= \mu\_{0} \frac{\partial \mu\_{0}}{\partial a}, \quad \mathcal{A}\_{1} = \mu\_{0} \frac{\partial \mu\_{1}}{\partial a} + \mu\_{1} \frac{\partial \mu\_{0}}{\partial a}, \quad \mathcal{B}\_{0} = \mu\_{0} \frac{\partial \nu\_{0}}{\partial \beta}, \quad \mathcal{B}\_{1} = \mu\_{0} \frac{\partial \nu\_{1}}{\partial \beta} + \mu\_{1} \frac{\partial \nu\_{0}}{\partial \beta}, \\ \mathcal{C}\_{0} &= \nu\_{0} \frac{\partial \mu\_{0}}{\partial a}, \quad \mathcal{C}\_{1} = \nu\_{0} \frac{\partial \mu\_{1}}{\partial a} + \nu\_{1} \frac{\partial \mu\_{0}}{\partial a}, \end{split}$$

Thus, we can easily obtain the recursive relationship by comparing two sides of Equation (33):

$$
\mu\_0(\alpha, \eta) = \frac{1}{2} - 8 \tanh(-2a), \quad \nu\_0(a, \eta) = 16 - 16 \tanh^2(-2a),
$$

For *m* = 0,

$$\mu\_1(\mathfrak{a}, \eta) = -8 \sec h^2 (-2a) \frac{\eta^\delta}{\Gamma(\delta + 1)}, \quad \nu\_1(\mathfrak{a}, \eta) = -32 \sec h^2 (-2a) \tanh(-2a) \frac{\eta^\delta}{\Gamma(\delta + 1)}.$$

For *m* = 1,

$$\begin{aligned} \mu\_2(a,\eta) &= -16 \sec h^2 (-2a) \left( 4 \sec h^2 (-2a) - 8 \tanh^2 (-2a) + 3 \tanh(-2a) \right) \frac{\eta^{2\delta}}{\Gamma(2\delta + 1)}, \\\nu\_2(a,\eta) &= -32 \sec h^2 (-2a) \{ 40 \sec h^2 (-2a) \tanh(-2a) + 96 \tanh(-2a) - 2 \tanh^2 (-2a) - 32 \tanh^3 (-2a) \}, \\\ &- 25 \sec h^2 (-2a) \{ \frac{\eta^{2\delta}}{\Gamma(2\delta + 1)}, \end{aligned}$$

In the same procedure, the remaining *μ<sup>m</sup>* and *ν<sup>m</sup>* (*m* ≥ 2) components of the NDM solution can be obtained smoothly. We therefore determine the sequence of alternatives as:

$$\mu(\mathfrak{a},\eta) = \sum\_{m=0}^{\infty} \mu\_m(\mathfrak{a},\mathfrak{\beta}) = \mu\_0(\mathfrak{a},\mathfrak{\beta}) + \mu\_1(\mathfrak{a},\mathfrak{\beta}) + \mu\_2(\mathfrak{a},\mathfrak{\beta}) + \mu\_3(\mathfrak{a},\mathfrak{\beta}) + \dots + \mu\_n(\mathfrak{a},\mathfrak{\beta})$$

$$\nu(\mathfrak{a},\eta) = \sum\_{m=0}^{\infty} \nu\_m(\mathfrak{a},\mathfrak{\beta}) = \nu\_0(\mathfrak{a},\mathfrak{\beta}) + \nu\_1(\mathfrak{a},\mathfrak{\beta}) + \nu\_2(\mathfrak{a},\mathfrak{\beta}) + \nu\_3(\mathfrak{a},\mathfrak{\beta}) + \dotsb$$

$$\begin{split} \mu(a,\eta) &= \frac{1}{2} - 8\tanh(-2a) - 8\sec h^2(-2a)\frac{\eta^\delta}{\Gamma(\delta+1)} \\ &- 16\sec h^2(-2a) \left( 4\sec h^2(-2a) - 8\tanh^2(-2a) + 3\tanh(-2a) \right) \frac{\eta^{2\delta}}{\Gamma(2\delta+1)} - \cdots \\ \nu(a,\eta) &= 16 - 16\tanh^2(-2a) - 32\sec h^2(-2a)\tanh(-2a) \frac{\eta^\delta}{\Gamma(\delta+1)} \\ &- 32\sec h^2(-2a) \{40\sec h^2(-2a)\tanh(-2a) + 96\tanh(-2a) - 2\tanh^2(-2a) - 32\tanh^3(-2a)\} \\ &- 25\sec h^2(-2a) \{\frac{\eta^{2\delta}}{\Gamma(2\delta+1)} - \cdots \end{split}$$

In Figures 1 and 2, the exact and natural decomposition method (NDM) solutions at an integer-order *δ* = 1 are represented for both *μ*(*α*, *η*) and *ν*(*α*, *η*) of Example 1. It is observed that NDM solutions are in good contact with the exact solution of the problems. In Figures 3 and 4, various fractional-order solutions of Example 1, at different fractional-orders, *δ* = 1, 0.8, 0.6, 0.4 and *η* = 1 are plotted. It is investigated that for Example 1, the fractional-order solutions are convergent to an integer-order solution for both *μ*(*α*, *η*) and *ν*(*α*, *η*).

**Figure 1.** Exact and NDM solution of *μ*(*α*, *η*) at *δ* = 1.

**Figure 2.** Exact and NDM solution of *ν*(*α*, *η*) at *δ* = 1.

**Figure 3.** Exact and NDM solution of *μ*(*α*, *η*) at different fractional order *δ* = 1, 0.8, 0.6, 0.4, and *η* = 1.

**Figure 4.** Exact and NDM solution of *ν*(*α*, *η*) at different fractional order *δ* = 1, 0.8, 0.6, 0.4, and *η* = 1.

#### *5.1. q-Homotopy Analysis Transform Method*

The Example 1 approximate solution with the help of **q-HATM**. After the Laplace transformation of Equation (29), we get:

$$\begin{split} \mathcal{E}\left\{\mu(a,\eta)\right\} &= \frac{1}{s} \left\{\mu(a,0)\right\} - \frac{1}{s^{\delta}} \mathcal{E}\left[\mu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} + \frac{\partial\mu(a,\eta)}{\partial a} + \frac{\partial\nu(a,\eta)}{\partial a}\right], \\ \mathcal{E}\left\{\nu(a,\eta)\right\} &= \frac{1}{s} \left\{\nu(a,0)\right\} - \frac{1}{s^{\delta}} \mathcal{E}\left[\mu(a,\eta)\frac{\partial\nu(a,\eta)}{\partial a} + \nu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} + 3\frac{\partial^{3}\mu(a,\eta)}{\partial a^{3}} - \frac{\partial^{2}\nu(a,\eta)}{\partial a^{2}}\right]. \end{split} \tag{34}$$

By the help of Equation (34) we define the nonlinear operator as:

$$\begin{split} &N^{\mathbb{L}}\left[\phi\_{1}(a,\eta;q),\phi\_{2}(a,\eta;q)\right] \\ &=\mathcal{L}\left[\phi\_{1}(a,\eta;q)-\frac{1}{s}\left\{\phi\_{1}(a,0)\right\}+\frac{1}{s^{\mathcal{S}}}\left\{\phi\_{1}(a,\eta)\frac{\partial\phi\_{1}(a,\eta)}{\partial a}+\frac{\partial\phi\_{1}(a,\eta)}{\partial a}+\frac{\partial\phi\_{2}(a,\eta)}{\partial a}\right\}\right], \\ &N^{\mathbb{L}}[\phi\_{1}(a,\eta;q),\phi\_{2}(a,\eta;q)] \\ &=\mathcal{L}\left[\phi\_{2}(a,\eta;q)-\frac{1}{s}\left\{\phi\_{2}(a,0)\right\}+\frac{1}{s^{\mathcal{S}}}\left\{\phi\_{1}(a,\eta)\frac{\partial\phi\_{2}(a,\eta)}{\partial a}+\phi\_{2}(a,\eta)\frac{\partial\phi\_{1}(a,\eta)}{\partial a}+3\frac{\partial^{3}\phi\_{1}(a,\eta)}{\partial a^{3}}-\frac{\partial^{2}\phi\_{2}(a,\eta)}{\partial a^{2}}\right\}\right]. \end{split} \tag{35}$$

By applying proposed algorithm, the deformation equation of *m*-th order is given as:

$$\begin{aligned} \mathbb{E}[\mu\_m(\boldsymbol{a}, \boldsymbol{\eta}) - \boldsymbol{K}\_m \boldsymbol{\mu}\_{m-1}(\boldsymbol{a}, \boldsymbol{\eta})] &= h \mathbb{R}\_{1,m}[\overrightarrow{\boldsymbol{\mu}}\_{m-1}, \overrightarrow{\boldsymbol{\nabla}}\_{m-1}], \\ \mathbb{E}[\boldsymbol{\nu}\_m(\boldsymbol{a}, \boldsymbol{\eta}) - \boldsymbol{K}\_m \boldsymbol{\nu}\_{m-1}(\boldsymbol{a}, \boldsymbol{\eta})] &= h \mathbb{R}\_{2,m}[\overrightarrow{\boldsymbol{\mu}}\_{m-1}, \overrightarrow{\boldsymbol{\nabla}}\_{m-1}], \end{aligned} \tag{36}$$

$$\begin{split} \vert \theta\_{1,m} \vert \overrightarrow{\boldsymbol{\mu}}\_{m-1}, \overrightarrow{\boldsymbol{\nu}}\_{m-1} \vert &= \mathcal{L} [\mu\_{m-1}(a, \eta) - (1 - \frac{k\_m}{n})\frac{1}{s} \{\frac{1}{2} - 8 \tanh(-2a)\} \\ &+ \frac{1}{s^2} \{\sum\_{j=0}^{m-1} \mu\_j(a, \eta) \frac{\partial \mu\_{m-1-j}(a, \eta)}{\partial a} + \frac{\partial \mu\_{m-1}(a, \eta)}{\partial a} + \frac{\partial \upsilon\_{m-1}(a, \eta)}{\partial a} \} \vert, \\ \vert \theta\_{2,m} \vert \overleftarrow{\boldsymbol{\mu}}\_{m-1}, \overrightarrow{\boldsymbol{\nu}}\_{m-1} \vert &= \mathcal{L} [\upsilon\_{m-1}(a, \eta) - \frac{1}{s} \{16 - 16 \tanh^2(-2a)\} \\ &+ \frac{1}{s^2} \{\sum\_{j=0}^{m-1} \mu\_j(a, \eta) \frac{\partial \upsilon\_{m-1-j}(a, \eta)}{\partial a} + \sum\_{j=0}^{m-1} \upsilon\_j(a, \eta) \frac{\partial \mu\_{m-j-1}(a, \eta)}{\partial a} + 3 \frac{\partial^3 \mu\_{m-1}(a, \eta)}{\partial a^3} - \frac{\partial^2 \upsilon\_{m-1}(a, \eta)}{\partial a^2} \} \vert. \end{split} \tag{37}$$

By applying inverse Laplace transform on Equation (36), we get:

$$\begin{aligned} \mu\_m(a,\eta) &= K\_m \mu\_{m-1}(a,\eta) + hL^{-1} \mathfrak{R}\_{1,m}[\overrightarrow{\mu}\_{m-1}, \overrightarrow{\nu}\_{m-1}],\\ \nu\_m(a,\eta) &= K\_m \nu\_{m-1}(a,\eta)] + hL^{-1} \mathfrak{R}\_{2,m}[\overrightarrow{\mu}\_{m-1}, \overrightarrow{\nu}\_{m-1}], \end{aligned} \tag{38}$$

By the help of given initial condition, we have:

$$\begin{aligned} \mu\_0(a,\eta) &= \frac{1}{2} - 8 \tanh(-2a), \\ \nu\_0(a,\eta) &= 16 - 16 \tanh^2(-2a). \end{aligned} \tag{39}$$

To find the value of *μ*0(*α*, *η*) and *ν*0(*α*, *η*), set *m* = 1 in Equation (38), then we get:

$$\begin{aligned} \mu\_1(\mathfrak{a}, \eta) &= K\_1 \mu\_0(\mathfrak{a}, \eta) + h \mathbb{E}^{-1} \mathfrak{R}\_{1,1} [\overrightarrow{\mu}^\flat \underset{\eta}{\cdot} \overrightarrow{\nu}^\flat ] , \\ \nu\_1(\mathfrak{a}, \eta) &= K\_1 \nu\_0(\mathfrak{a}, \eta) [+h \mathbb{E}^{-1} \mathfrak{R}\_{2,1} [\overrightarrow{\mu}^\flat \underset{\eta}{\cdot} \overrightarrow{\nu}^\flat ] , \end{aligned} \tag{40}$$

From Equation (37) for *m* = 1, we get:

$$\begin{split} \mathfrak{R}\_{1,1}[\overrightarrow{\mu}\_{0},\overrightarrow{\nu}\_{0}] &= \mathfrak{E}[\mu\_{0}(a,\eta)] - (1 - \frac{k\_{1}}{n})\frac{1}{s} \{\frac{1}{2} - 8\tanh(-2a)\} \\ &+ \frac{1}{s^{3}}\mathfrak{E}[\{\mu\_{0}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha} + \frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha} + \frac{\partial\nu\_{0}(a,\eta)}{\partial\alpha}\}], \\ \mathfrak{R}\_{2,1}[\overrightarrow{\mu}\_{0},\overrightarrow{\nu}\_{0}] &= \mathfrak{E}[\nu\_{0}(a,\eta)] - (1 - \frac{k\_{1}}{n})\frac{1}{s} \{16 - 16\tanh^{2}(-2a)\} \\ &+ \frac{1}{s^{3}}\mathfrak{E}[\{\mu\_{0}(a,\eta)\frac{\partial\nu\_{0}(a,\eta)}{\partial\alpha} + \nu\_{0}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha} + 3\frac{\partial^{3}\mu\_{0}(a,\eta)}{\partial\alpha^{3}} - \frac{\partial^{2}\nu\_{0}(a,\eta)}{\partial\alpha^{2}}\}]. \end{split} \tag{41}$$

Then by using Equations (25) and (41) in Equation (40), we get:

$$\begin{split} \mu\_{1}(a,\eta) &= h\mathcal{L}^{-1}[\frac{1}{s}\{\frac{1}{2} - 8\tanh(-2a)\} - (1 - \frac{0}{n})\frac{1}{s}\{\frac{1}{2} - 8\tanh(-2a)\} \\ &+ \frac{1}{s^{\delta}}\mathcal{L}[\{\mu\_{0}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha} + \frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha} + \frac{\partial\nu\_{0}(a,\eta)}{\partial\alpha}\}]), \\ \nu\_{1}(a,\eta) &= h\mathcal{L}^{-1}[\frac{1}{s}\{16 - 16\tanh^{2}(-2a)\} - (1 - \frac{0}{n})\frac{1}{s}\{16 - 16\tanh^{2}(-2a)\} \\ &+ \frac{1}{s^{\delta}}\mathcal{L}[\{\mu\_{0}(a,\eta)\frac{\partial\nu\_{0}(a,\eta)}{\partial\alpha} + \nu\_{0}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha} + 3\frac{\partial^{3}\mu\_{0}(a,\eta)}{\partial a^{3}} - \frac{\partial^{2}\nu\_{0}(a,\eta)}{\partial a^{2}}\}]), \\ \mu\_{1}(a,\eta) &= -8h\operatorname{sech}^{2}(-2a)\frac{\eta^{\delta}}{\Gamma(\delta+1)}, \quad \nu\_{1}(a,\eta) = -32h\operatorname{sech}^{2}(-2a)\tanh(-2a)\frac{\eta^{\delta}}{\Gamma(\delta+1)}. \end{split} \tag{10.51}$$

Similarly from Equations (40) and (41) for *m* = 2, we have:

$$\begin{split} \mu\_{2}(a,\eta) &= n\mu\_{1}(a,\eta) + \hbar\ell^{-1}[\mathcal{E}[\mu\_{1}(a,\eta)] - (1 - \frac{n}{n})\frac{1}{s}(\frac{1}{2} - 8\tanh(-2a)) + \frac{1}{s^{2}}\mathcal{E}[\{\mu\_{0}(a,\eta)\frac{\partial\mu\_{1}(a,\eta)}{\partial\mathfrak{a}}\}] \\ &+ \mu\_{1}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\mathfrak{a}} + \frac{\partial\mu\_{1}(a,\eta)}{\partial\mathfrak{a}} + \frac{\partial\nu\_{1}(a,\eta)}{\partial\mathfrak{a}})\}], \\ \nu\_{2}(a,\eta) &= m\nu\_{1}(a,\eta) + \hbar\mathcal{L}^{-1}[\mathcal{E}[\mu\_{1}(a,\eta)] - (1 - \frac{n}{n})\frac{1}{s}(16 - 16\tanh^{2}(-2a)) + \frac{1}{s^{2}}\mathcal{E}[\{\mu\_{0}(a,\eta)\frac{\partial\mu\_{1}(a,\eta)}{\partial\mathfrak{a}}\}] \\ &+ \mu\_{1}(a,\eta)\frac{\partial\nu\_{0}(a,\eta)}{\partial\mathfrak{a}} + \nu\_{0}(a,\eta)\frac{\partial\mu\_{1}(a,\eta)}{\partial\mathfrak{a}} + \nu\_{1}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\mathfrak{a}} + 3\frac{3^{3}\mu\_{1}(a,\eta)}{\partial\mathfrak{a}^{3}} - \frac{\partial^{2}\nu\_{1}(a,\eta)}{\partial\mathfrak{a}^{2}}\})]. \end{split} \tag{43}$$

In the case of simplified, the above calculation eliminates as described:

$$\begin{split} \mu\_{2}(\boldsymbol{a},\boldsymbol{\eta}) &= -8(\boldsymbol{n}+\boldsymbol{h})\hbar \sec \boldsymbol{h}^{2}(-2\boldsymbol{a}) \frac{\eta^{\delta}}{\Gamma(\delta+1)} - 16\boldsymbol{h}^{2} \sec \boldsymbol{h}^{2}(-2\boldsymbol{a})(4\sec \boldsymbol{h}^{2}(-2\boldsymbol{a}) \\ &- 8\tanh^{2}(-2\boldsymbol{a}) + 3\tanh(-2\boldsymbol{a})) \frac{\eta^{2\delta}}{\Gamma(2\delta+1)}, \\ \nu\_{2}(\boldsymbol{a},\boldsymbol{\eta}) &= -32(\boldsymbol{n}+\boldsymbol{h})\hbar \sec \boldsymbol{h}^{2}(-2\boldsymbol{a})\tanh(-2\boldsymbol{a})\frac{\eta^{\delta}}{\Gamma(\delta+1)} - 32\boldsymbol{h}^{2} \sec \boldsymbol{h}^{2}(-2\boldsymbol{a})\{40\sec \boldsymbol{h}^{2}(-2\boldsymbol{a})\tanh(-2\boldsymbol{a})\} \\ &+ 96\tanh(-2\boldsymbol{a}) - 2\tanh^{2}(-2\boldsymbol{a}) - 32\tanh^{3}(-2\boldsymbol{a}) - 25\sec \boldsymbol{h}^{2}(-2\boldsymbol{a})\left\{\frac{\eta^{2\delta}}{\Gamma(2\delta+1)}\right\} \end{split}$$

The rest of the iterative terms can be used in the same way. Formerly, the family of q-homotopy analysis transform technique series result of Equation (29) is assumed by:

$$\begin{aligned} \mu(\mathfrak{a}, \eta) &= \mu\_0(\mathfrak{a}, \eta) + \sum\_{m=1}^{\infty} \mu\_m(\mathfrak{a}, \eta) (\frac{1}{n})^m, \\ \nu(\mathfrak{a}, \eta) &= \nu\_0(\mathfrak{a}, \eta) + \sum\_{m=1}^{\infty} \nu\_m(\mathfrak{a}, \eta) (\frac{1}{n})^m, \end{aligned} \tag{44}$$

The exact solution of Equation (29) at *δ* = 1,

$$\begin{aligned} \mu(\mathfrak{a}, \eta) &= \frac{1}{2} - 8 \tanh \left\{ -2 \left( \mathfrak{a} - \frac{\eta}{2} \right) \right\}, \\ \nu(\mathfrak{a}, \eta) &= 16 - 16 \tanh^2 \left\{ -2 \left( \mathfrak{a} - \frac{\eta}{2} \right) \right\}. \end{aligned} \tag{45}$$

In Figure 5, the graph of exact and q-HATM solutions for *μ*(*α*, *η*) of Example 1 are displayed. It is observed that, the solutions of q-HATM are in good agreement with the exact and NDM solutions. Similarly Figure 6, express the exact and q-HATM solutions for *ν*(*α*, *η*). The plot representation also confirmed the higher accuracy of the proposed method with the exact solution for *ν*(*α*, *η*). Furthermore, the graphical representations of the solutions of the proposed method have reflected its applicability and reliability. This provides the motivation to apply the current techniques for other fractional-order partial differential equations.

**Figure 5.** Exact and q-HATM solution of *μ*(*α*, *η*) at *δ* = 1.

**Figure 6.** Exact and q-HATM solution of *ν*(*α*, *η*) at *δ* = 1.

$$\begin{split} &D^{\delta}\_{\eta}\mu(a,\eta) + \mu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} + \frac{1}{2}\frac{\partial\mu(a,\eta)}{\partial a} + \frac{\partial\nu(a,\eta)}{\partial a} = 0, \\ &D^{\delta}\_{\eta}\nu(a,\eta) + \mu(a,\eta)\frac{\partial\nu(a,\eta)}{\partial a} + \nu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} - \frac{1}{2}\frac{\partial^{2}\nu(a,\eta)}{\partial a^{2}} = 0, \\ &0 < \delta \le 1, \quad 0 < \eta \le 1, \quad -100 \le a \le 100, \end{split} \tag{46}$$

*with the initial condition:*

$$\begin{cases} & \mu(\mathfrak{a},0) = \mathfrak{z} - \mathfrak{x}\coth[\mathfrak{x}(\mathfrak{a}+\theta)], \\ & \nu(\mathfrak{a},0) = -\mathfrak{x}^2\mathrm{cosech}^2[\mathfrak{x}(\mathfrak{a}+\theta)]. \end{cases} \tag{47}$$

*Firstly, we will solve this scheme by using the NDM.*

After the natural transformation of Equation (46), we get:

*N*<sup>+</sup> *∂δμ*(*α*, *η*) *∂η<sup>δ</sup>* <sup>=</sup> <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂μ*(*α*, *η*) *∂α* <sup>+</sup> 1 2 *∂μ*(*α*, *η*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* " , *N*<sup>+</sup> *∂δν*(*α*, *η*) *∂η<sup>δ</sup>* <sup>=</sup> <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂ν*(*α*, *η*) *∂α* <sup>+</sup> *<sup>ν</sup>*(*α*, *<sup>η</sup>*) *∂μ*(*α*, *η*) *∂α* <sup>−</sup> <sup>1</sup> 2 *∂*2*ν*(*α*, *η*) *∂α*<sup>2</sup> " , *sδ <sup>u</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> {*μ*(*α*, *<sup>η</sup>*)} <sup>−</sup> *<sup>s</sup>δ*−<sup>1</sup> *<sup>u</sup><sup>δ</sup> <sup>μ</sup>*(*α*, 0) = <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂μ*(*α*, *η*) *∂α* <sup>+</sup> 1 2 *∂μ*(*α*, *η*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* " *sδ <sup>u</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> {*ν*(*α*, *<sup>η</sup>*)} <sup>−</sup> *<sup>s</sup>δ*−<sup>1</sup> *<sup>u</sup><sup>δ</sup> <sup>ν</sup>*(*α*, 0) = <sup>−</sup>*N*<sup>+</sup> ! *μ*(*α*, *η*) *∂ν*(*α*, *η*) *∂α* <sup>+</sup> *<sup>ν</sup>*(*α*, *<sup>η</sup>*) *∂μ*(*α*, *η*) *∂α* <sup>−</sup> <sup>1</sup> 2 *∂*2*ν*(*α*, *η*) *∂α*<sup>2</sup> " ,

The above algorithm is reduced to the simplified form as:

$$\begin{split} \{N^{+}\left\{\mu(a,\eta)\right\}\} &= \frac{1}{\mathcal{S}}\left\{\mu(a,0)\right\} - \frac{\mathsf{u}^{\delta}}{\mathsf{s}^{\delta}}\mathcal{N}^{+}\left[\mu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} + \frac{1}{2}\frac{\partial\mu(a,\eta)}{\partial a} + \frac{\partial\nu(a,\eta)}{\partial a}\right],\\ \{N^{+}\left\{\nu(a,\eta)\right\}\} &= \frac{1}{\mathcal{S}}\left\{\nu(a,0)\right\} - \frac{\mathsf{u}^{\delta}}{\mathsf{s}^{\delta}}\mathcal{N}^{+}\left[\mu(a,\eta)\frac{\partial\nu(a,\eta)}{\partial a} + \nu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} - \frac{1}{2}\frac{\partial^{2}\nu(a,\eta)}{\partial a^{2}}\right]. \end{split} \tag{48}$$

Applying the inverse natural transformation, we get:

$$\begin{split} \mu(\mathfrak{a},\eta) &= \mu(\mathfrak{a},0) - N^- \left[ \frac{\mathfrak{a}^\delta}{s^\delta} N^+ \left[ \mu(\mathfrak{a},\eta) \frac{\partial \mu(\mathfrak{a},\eta)}{\partial \mathfrak{a}} + \frac{1}{2} \frac{\partial \mu(\mathfrak{a},\eta)}{\partial \mathfrak{a}} + \frac{\partial \nu(\mathfrak{a},\eta)}{\partial \mathfrak{a}} \right] \right], \\ \nu(\mathfrak{a},\eta) &= \nu(\mathfrak{a},0) - N^- \left[ \frac{\mathfrak{a}^\delta}{s^\delta} N^+ \left[ \mu(\mathfrak{a},\eta) \frac{\partial \nu(\mathfrak{a},\eta)}{\partial \mathfrak{a}} + \nu(\mathfrak{a},\eta) \frac{\partial \mu(\mathfrak{a},\eta)}{\partial \mathfrak{a}} - \frac{1}{2} \frac{\partial^2 \nu(\mathfrak{a},\eta)}{\partial \mathfrak{a}^2} \right] \right], \end{split} \tag{49}$$

Assume that the unknown functions *μ*(*α*, *η*) and *ν*(*α*, *η*) infinite series solution is as follows:

$$
\mu(\mathfrak{a}, \eta) = \sum\_{m=0}^{\infty} \mu\_m(\mathfrak{a}, \eta), \quad \text{and} \quad \nu(\mathfrak{a}, \eta) = \sum\_{m=0}^{\infty} \nu\_m(\mathfrak{a}, \eta),
$$

Remember that *μμα* = ∑<sup>∞</sup> *<sup>m</sup>*=<sup>0</sup> <sup>A</sup>*m*, *μνα* <sup>=</sup> <sup>∑</sup><sup>∞</sup> *<sup>m</sup>*=<sup>0</sup> <sup>B</sup>*<sup>m</sup>* and *νμα* <sup>=</sup> <sup>∑</sup><sup>∞</sup> *<sup>m</sup>*=<sup>0</sup> C*<sup>m</sup>* are the Adomian polynomials and the nonlinear terms were characterized. Using such terms, Equation (49) can be rewritten in the form:

∞ ∑ *m*=0 *μm*(*α*, *η*) = *μ*(*α*, 0) − *N*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 A*<sup>m</sup>* + 1 2 *∂μ*(*α*, *η*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* , ∞ ∑ *m*=0 *νm*(*α*, *η*) = *ν*(*α*, 0) − *N*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 B*<sup>m</sup>* + ∞ ∑ *m*=0 <sup>C</sup>*<sup>m</sup>* <sup>−</sup> <sup>1</sup> 2 *∂*2*ν*(*α*, *η*) *∂α*<sup>2</sup> , ∞ ∑ *m*=0 *μm*(*α*, *η*) = *ξ* − *κ* coth[*κ*(*α* + *θ*)] − *N*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 A*<sup>m</sup>* + 1 2 *∂μ*(*α*, *η*) *∂α* <sup>+</sup> *∂ν*(*α*, *η*) *∂α* , ∞ ∑ *m*=0 *<sup>ν</sup>m*(*α*, *<sup>η</sup>*) = <sup>−</sup>*κ*2*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] <sup>−</sup> *<sup>N</sup>*<sup>−</sup> *uδ <sup>s</sup><sup>δ</sup> <sup>N</sup>*<sup>+</sup> ∞ ∑ *m*=0 B*<sup>m</sup>* + ∞ ∑ *m*=0 <sup>C</sup>*<sup>m</sup>* <sup>−</sup> <sup>1</sup> 2 *∂*2*ν*(*α*, *η*) *∂α*<sup>2</sup> , (50)

According to Equation (7), all forms of non-linearity the Adomian polynomials can be defined as:

$$\begin{split} \mathcal{A}\_{0} &= \mu\_{0} \frac{\partial \mu\_{0}}{\partial \alpha}, \quad \mathcal{A}\_{1} = \mu\_{0} \frac{\partial \mu\_{1}}{\partial \alpha} + \mu\_{1} \frac{\partial \mu\_{0}}{\partial \alpha}, \quad \mathcal{B}\_{0} = \mu\_{0} \frac{\partial \upsilon\_{0}}{\partial \beta}, \quad \mathcal{B}\_{1} = \mu\_{0} \frac{\partial \upsilon\_{1}}{\partial \beta} + \mu\_{1} \frac{\partial \upsilon\_{0}}{\partial \beta}, \\ \mathcal{C}\_{0} &= \upsilon\_{0} \frac{\partial \mu\_{0}}{\partial \alpha'}, \quad \mathcal{C}\_{1} = \upsilon\_{0} \frac{\partial \mu\_{1}}{\partial \alpha} + \upsilon\_{1} \frac{\partial \mu\_{0}}{\partial \alpha'}, \end{split}$$

Thus, we can easily obtain the recursive relationship by comparing two sides of Equation (50):

$$
\mu\_0(\mathfrak{a}, \eta) = \mathfrak{z} - \kappa \coth[\kappa(\mathfrak{a} + \theta)], \quad \nu\_0(\mathfrak{a}, \eta) = -\kappa^2 \alpha \mathrm{sech}^2[\kappa(\mathfrak{a} + \theta)],
$$

For *m* = 0,

$$\begin{aligned} \mu\_1(\mathfrak{a}, \eta) &= -\lg^2 \text{coscech}^2[\kappa(\mathfrak{a} + \theta)] \frac{\eta^\delta}{\Gamma(\delta + 1)}, \\ \nu\_1(\mathfrak{a}, \eta) &= -\lg^2 \text{coscech}^2[\kappa(\mathfrak{a} + \theta)] \coth[\kappa(\mathfrak{a} + \theta)] \frac{\eta^\delta}{\Gamma(\delta + 1)}. \end{aligned}$$

For *m* = 1,

$$\begin{split} \mu\_{2}(\boldsymbol{a},\boldsymbol{\eta}) &= \xi \kappa^{4} \mathrm{cscech}^{2}[\kappa(\boldsymbol{a}+\boldsymbol{\theta})] \left\{ \frac{2\xi^{3}\kappa\Gamma(2\delta+1)\eta^{3\delta}}{(\Gamma(\delta+1))^{2}\Gamma(3\delta+1)} - \frac{(3\coth^{2}([\kappa(\boldsymbol{a}+\boldsymbol{\theta})]-1))\eta^{2\delta}}{\Gamma(2\delta+1)} \right\}, \\ \nu\_{2}(\boldsymbol{a},\boldsymbol{\eta}) &= \frac{1}{\Gamma(\delta+1)} [2\xi\kappa^{5}\mathrm{cscech}^{2}[\kappa(\boldsymbol{a}+\boldsymbol{\theta})]] [\frac{\xi\kappa\mathrm{cscech}^{2}(3\coth^{2}([\kappa(\boldsymbol{a}+\boldsymbol{\theta})]-1))\eta^{3\delta}}{\Gamma(\delta+1)\Gamma(3\delta+1)} \\ &+ \frac{2\xi\kappa\mathrm{cscech}^{2}\coth^{2}([\kappa(\boldsymbol{a}+\boldsymbol{\theta})])\eta^{3\delta}}{\Gamma(\delta+1)\Gamma(3\delta+1)} - \frac{2\xi\kappa\mathrm{coth}(3\csc\kappa\mathrm{h}^{2}([\kappa(\boldsymbol{a}+\boldsymbol{\theta})]-1))\eta^{2\delta}}{\Gamma(2\delta+1)}. \end{split}$$

In the same procedure, the remaining *μ<sup>m</sup>* and *ν<sup>m</sup>* (*m* ≥ 2) components of the NDM solution can be obtained smoothly. Thus, we determine the sequence of alternatives as:

*μ*(*α*, *η*) = ∞ ∑ *m*=0 *μm*(*α*, *β*) = *μ*0(*α*, *β*) + *μ*1(*α*, *β*) + *μ*2(*α*, *β*) + *μ*3(*α*, *β*) + ··· *ν*(*α*, *η*) = ∞ ∑ *m*=0 *νm*(*α*, *β*) = *ν*0(*α*, *β*) + *ν*1(*α*, *β*) + *ν*2(*α*, *β*) + *ν*3(*α*, *β*) + ··· *<sup>μ</sup>*(*α*, *<sup>η</sup>*) = *<sup>ξ</sup>* <sup>−</sup> *<sup>κ</sup>* coth[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] <sup>−</sup> *ξκ*2*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] *<sup>η</sup><sup>δ</sup>* Γ(*δ* + 1) <sup>+</sup> *ξκ*4*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] - 2*ξκ*Γ(2*δ* + 1)*η*3*<sup>δ</sup>* (Γ(*<sup>δ</sup>* <sup>+</sup> <sup>1</sup>))2Γ(3*<sup>δ</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> (3 coth2([*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] <sup>−</sup> <sup>1</sup>))*η*2*<sup>δ</sup>* Γ(2*δ* + 1) −··· *<sup>ν</sup>*(*α*, *<sup>η</sup>*) = <sup>−</sup>*κ*2*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] <sup>−</sup> *ξκ*2*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] coth[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] *<sup>η</sup><sup>δ</sup>* Γ(*δ* + 1)

$$\begin{split} &+\frac{1}{\Gamma(\delta+1)}[2\xi\mathsf{x}^{5}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{h}^{2}(\mathsf{x}(a+\theta))] \frac{\xi\mathsf{x}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{h}^{2}(3\mathsf{c}\mathsf{c}\mathsf{h}^{2}([\mathsf{x}(a+\theta)]-1))\eta^{3\delta}}{\Gamma(\delta+1)\Gamma(3\delta+1)} \\ &+\frac{2\xi\mathsf{x}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{h}^{2}\mathsf{c}\mathsf{c}\mathsf{h}^{2}([\mathsf{x}(a+\theta)])\eta^{3\delta}}{\Gamma(\delta+1)\Gamma(3\delta+1)} - \frac{2\xi\mathsf{x}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{h}(3\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{c}\mathsf{h}^{2}([\mathsf{x}(a+\theta)]-1))\eta^{2\delta}}{\Gamma(2\delta+1)} - \cdots \end{split}$$

Figures 7 and 8 describe the graphical behavior of both the unknown variables *μ*(*α*, *η*) and *ν*(*α*, *η*) of Example 2 at an integer-order *δ* = 1 respectively. The procedures of NDM and q-HATM are implemented to obtain the desire accuracy. The higher accuracy and rate of convergence are achieved by the proposed techniques as shown in Figure 9. The plot analysis demonstrates the validity and accuracy of the proposed techniques and considered to be the best techniques to solve other fractional-order problems.

**Figure 7.** Exact and NDM solution of *μ*(*α*, *η*) at *δ* = 1.

**Figure 8.** Exact and NDM solution of *ν*(*α*, *η*) at *δ* = 1.

**Figure 9.** Error plot of *μ*(*α*, *η*) and *ν*(*α*, *η*).

#### *5.2. q-Homotopy Analysis Transform Method*

The Example 1 approximate solution with the help of **q-HATM**. By taking the Laplace transformation of Equation (46), we get:

$$\begin{split} \mathcal{L}\left\{\mu(a,\eta)\right\} &= \frac{1}{s} \left\{\mu(a,0)\right\} - \frac{1}{s^{\delta}} \mathcal{E}\left[\mu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} + \frac{1}{2}\frac{\partial\mu(a,\eta)}{\partial a} + \frac{\partial\nu(a,\eta)}{\partial a}\right], \\ \mathcal{E}\left\{\nu(a,\eta)\right\} &= \frac{1}{s} \left\{\nu(a,0)\right\} - \frac{1}{s^{\delta}} \mathcal{E}\left[\mu(a,\eta)\frac{\partial\nu(a,\eta)}{\partial a} + \nu(a,\eta)\frac{\partial\mu(a,\eta)}{\partial a} - \frac{1}{2}\frac{\partial^{2}\nu(a,\eta)}{\partial a^{2}}\right], \end{split} \tag{51}$$

Using Equation (51) we define the nonlinear operator as:

$$\begin{split} &N^{1}\left[\phi\_{1}(a,\eta;q),\phi\_{2}(a,\eta;q)\right] \\ &=\mathcal{L}\left[\phi\_{1}(a,\eta;q)-\frac{1}{s}\left\{\phi\_{1}(a,0)\right\}+\frac{1}{s^{\delta}}\left\{\phi\_{1}(a,\eta)\frac{\partial\phi\_{1}(a,\eta)}{\partial a}+\frac{1}{2}\frac{\partial\phi\_{1}(a,\eta)}{\partial a}+\frac{\partial\phi\_{2}(a,\eta)}{\partial a}\right\}\right], \\ &N^{2}[\phi\_{1}(a,\eta;q),\phi\_{2}(a,\eta;q)] \\ &=\mathcal{L}\left[\phi\_{2}(a,\eta;q)-\frac{1}{s}\left\{\phi\_{2}(a,0)\right\}+\frac{1}{s^{\delta}}\left\{\phi\_{1}(a,\eta)\frac{\partial\phi\_{2}(a,\eta)}{\partial a}+\phi\_{2}(a,\eta)\frac{\partial\phi\_{1}(a,\eta)}{\partial a}-\frac{1}{2}\frac{\partial^{2}\phi\_{2}(a,\eta)}{\partial a^{2}}\right\}\right]. \end{split} \tag{52}$$

By applying the proposed algorithm, the deformation equation of *m*-th order is given as:

$$\begin{aligned} \mathbb{E}[\mu\_m(\mathfrak{a}, \eta) - \mathsf{K}\_m \mathsf{u}\_{m-1}(\mathfrak{a}, \eta)] &= h \mathfrak{R}\_{1,m}[\overrightarrow{\mu}\_{m-1}, \overrightarrow{\nabla}\_{m-1}]\_{\prime} \\ \mathbb{E}[\upsilon\_m(\mathfrak{a}, \eta) - \mathsf{K}\_m \mathsf{v}\_{m-1}(\mathfrak{a}, \eta)] &= h \mathfrak{R}\_{2,m}[\overrightarrow{\mu}\_{m-1}, \overrightarrow{\nabla}\_{m-1}]\_{\prime} \end{aligned} \tag{53}$$

$$\begin{split} \mathfrak{R}\_{1,m}[\overrightarrow{\mu}\_{m-1}, \overrightarrow{\nu}\_{m-1}] &= \mathfrak{L}[\mu\_{m-1}(a, \eta) - (1 - \frac{k\_m}{n})\frac{1}{s} \{\xi - \kappa \coth[\kappa(a + \theta)]\} \\ &+ \frac{1}{s^\delta} \{\sum\_{j=0}^{m-1} \mu\_j(a, \eta) \frac{\partial \mu\_{m-1-j}(a, \eta)}{\partial \alpha} + \frac{1}{2} \frac{\partial \mu\_{m-1}(a, \eta)}{\partial \alpha} + \frac{\partial \nu\_{m-1}(a, \eta)}{\partial \alpha} \}), \\ \mathfrak{R}\_{2,m}[\overleftarrow{\mu}\_{m-1}, \overrightarrow{\nu}\_{m-1}] &= \mathfrak{L}[\nu\_{m-1}(a, \eta) - \frac{1}{s} \{-\kappa^2 \text{csech}^2[\kappa(a + \theta)]\} \\ &+ \frac{1}{s^\delta} \{\sum\_{j=0}^{m-1} \mu\_j(a, \eta) \frac{\partial \nu\_{m-1-j}(a, \eta)}{\partial \alpha} + \sum\_{j=0}^{m-1} \nu\_j(a, \eta) \frac{\partial \mu\_{m-j-1}(a, \eta)}{\partial \alpha} - \frac{1}{2} \frac{\partial^2 \nu\_{m-1}(a, \eta)}{\partial \alpha^2} \}). \end{split} \tag{54}$$

By applying the inverse Laplace transform on Equation (53), we get:

$$\begin{aligned} \mu\_{\mathfrak{m}}(a,\eta) &= K\_{\mathfrak{m}}\mu\_{\mathfrak{m}-1}(a,\eta) + hL^{-1}\mathfrak{R}\_{1,\mathfrak{m}}[\overrightarrow{\mu}\_{\mathfrak{m}-1}, \overrightarrow{\nu}\_{\mathfrak{m}-1}],\\ \nu\_{\mathfrak{m}}(a,\eta) &= K\_{\mathfrak{m}}\nu\_{\mathfrak{m}-1}(a,\eta)] + hL^{-1}\mathfrak{R}\_{2,\mathfrak{m}}[\overrightarrow{\mu}\_{\mathfrak{m}-1}, \overrightarrow{\nu}\_{\mathfrak{m}-1}].\end{aligned} \tag{55}$$

By the help of the given initial condition, we have:

$$\begin{aligned} \mu\_0(\boldsymbol{\alpha}, \boldsymbol{\eta}) &= \boldsymbol{\xi} - \kappa \coth[\kappa(\boldsymbol{\alpha} + \boldsymbol{\theta})], \\ \nu\_0(\boldsymbol{\alpha}, \boldsymbol{\eta}) &= -\kappa^2 \cosh^2[\kappa(\boldsymbol{\alpha} + \boldsymbol{\theta})]. \end{aligned} \tag{56}$$

To find the value of *μ*0(*α*, *η*) and *ν*0(*α*, *η*), set *m* = 1 in Equation (38), then we get:

$$\begin{aligned} \mu\_1(\mathfrak{a}, \eta) &= K\_1 \mu\_0(\mathfrak{a}, \eta) + h \mathbb{E}^{-1} \mathfrak{R}\_{1,1} [\overrightarrow{\mu}^\flat \llcorner \overrightarrow{\nu}^\flat \llcorner], \\ \nu\_1(\mathfrak{a}, \eta) &= K\_1 \nu\_0(\mathfrak{a}, \eta) [+h \mathbb{E}^{-1} \mathfrak{R}\_{2,1} [\overrightarrow{\mu}^\flat \llcorner \overrightarrow{\nu}^\flat \ll], \end{aligned} \tag{57}$$

From Equation (54) for *m* = 1, we conclude:

$$\begin{split} \mathfrak{R}\_{1,1}[\overrightarrow{\mu}\_{0}^{\dagger}\overrightarrow{\nu}\_{0}^{\prime}] &= \mathfrak{E}[\mu\_{0}(\mathsf{a},\mathsf{\eta})] - (1 - \frac{k\_{1}}{n})\frac{1}{s} \{\mathbb{S} - \kappa \coth[\kappa(\mathsf{a} + \theta)]\} \\ &+ \frac{1}{s^{\sharp}}\mathfrak{E}[\{\mu\_{0}(\mathsf{a},\mathsf{\eta})\frac{\partial \mu\_{0}(\mathsf{a},\mathsf{\eta})}{\partial \mathsf{a}} + \frac{1}{2}\frac{\partial \mu\_{0}(\mathsf{a},\mathsf{\eta})}{\partial \mathsf{a}} + \frac{\partial \upsilon\_{0}(\mathsf{a},\mathsf{\eta})}{\partial \mathsf{a}}\}], \\ \mathfrak{R}\_{2,1}[\overrightarrow{\mu}\_{0}^{\dagger}\overrightarrow{\nu}\_{0}] &= \mathfrak{E}[\upsilon\_{0}(\mathsf{a},\mathsf{\eta})] - (1 - \frac{k\_{1}}{n})\frac{1}{s} \{-\kappa^{2} \text{oscech}^{2}[\kappa(\mathsf{a} + \theta)]\} \\ &+ \frac{1}{s^{\sharp}}\mathfrak{E}[\{\mu\_{0}(\mathsf{a},\mathsf{\eta})\frac{\partial \upsilon\_{0}(\mathsf{a},\mathsf{\eta})}{\partial \mathsf{a}} + \upsilon\_{0}(\mathsf{a},\mathsf{\eta})\frac{\partial \mu\_{0}(\mathsf{a},\mathsf{\eta})}{\partial \mathsf{a}} - \frac{1}{2}\frac{\partial^{2} \upsilon\_{0}(\mathsf{a},\mathsf{\eta})}{\partial \mathsf{a}^{2}}\}, \end{split} \tag{58}$$

Then by using Equations (25) and (58) in Equation (57), we get

$$\begin{split} \mu\_{1}(a,\eta) &= h\mathcal{L}^{-1}[\frac{1}{s}\{\overline{\xi}-\mathsf{x}\coth[\mathsf{x}(a+\theta)]\}-(1-\frac{0}{n})\frac{1}{s}\{\overline{\xi}-\mathsf{x}\coth[\mathsf{x}(a+\theta)]\} \\ &+\frac{1}{s^{\mathcal{S}}}\mathcal{L}[\{\mu\_{0}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha}+\frac{1}{2}\frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha}+\frac{\partial\upsilon\_{0}(a,\eta)}{\partial\alpha}\}]), \\ \nu\_{1}(a,\eta) &= h\mathcal{L}^{-1}[\frac{1}{s}\{-\kappa^{2}\mathrm{cscech}^{2}[\kappa(a+\theta)]\}-(1-\frac{0}{n})\frac{1}{s}\{-\kappa^{2}\mathrm{cscech}^{2}[\kappa(a+\theta)]\} \\ &+\frac{1}{s^{\mathcal{S}}}\mathcal{L}[\{\mu\_{0}(a,\eta)\frac{\partial\upsilon\_{0}(a,\eta)}{\partial\alpha}+\upsilon\_{0}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\alpha}-\frac{1}{2}\frac{\partial^{2}\upsilon\_{0}(a,\eta)}{\partial\alpha^{2}}\}]), \\ \nu\_{1}(a,\eta) &= -\xi\hbar^{2}\mathrm{cscech}^{2}[\kappa(a+\theta)]\frac{\eta^{\delta}}{\Gamma(\delta+1)}, \quad \nu\_{1}(a,\eta) = -\xi\hbar\mathrm{cscech}^{2}[\kappa(a+\theta)]\coth[\kappa(a+\theta)]\frac{\eta^{\delta}}{\Gamma(\delta+1)}, \end{split} \tag{59}$$

Similarly from Equations (57) and (58) for *m* = 2, we have:

$$\begin{split} \mu\_{2}(a,\eta) &= n\mu\_{1}(a,\eta) + \hbar\mathcal{L}^{-1}[\mathcal{L}[\mu\_{1}(a,\eta)] - (1 - \frac{n}{n})\frac{1}{s} \{\sharp - \mathbf{x}\coth[\mathbf{x}(a+\theta)]\} + \frac{1}{s^{2}}\mathcal{L}[\{\mu\_{0}(a,\eta)\}\frac{\partial\mu\_{1}(a,\eta)}{\partial\mathbf{a}}] \\ &+ \mu\_{1}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\mathbf{a}} + \frac{1}{2}\frac{\partial\mu\_{1}(a,\eta)}{\partial\mathbf{a}} + \frac{\partial\nu\_{1}(a,\eta)}{\partial\mathbf{a}} \}) \}, \\ \nu\_{2}(a,\eta) &= n\nu\_{1}(a,\eta) + \hbar\mathcal{L}^{-1}[\mathcal{L}[\mu\_{1}(a,\eta)] - (1 - \frac{n}{n})\frac{1}{s} \{-\kappa^{2}\text{cscech}^{2}[\kappa(a+\theta)]\} + \frac{1}{s^{2}}\mathcal{L}[\{\mu\_{0}(a,\eta)\}\frac{\partial\nu\_{1}(a,\eta)}{\partial\mathbf{a}} \\ &+ \mu\_{1}(a,\eta)\frac{\partial\nu\_{0}(a,\eta)}{\partial\mathbf{a}} + \nu\_{0}(a,\eta)\frac{\partial\mu\_{1}(a,\eta)}{\partial\mathbf{a}} + \nu\_{1}(a,\eta)\frac{\partial\mu\_{0}(a,\eta)}{\partial\mathbf{a}} - \frac{1}{2}\frac{\partial^{2}\nu\_{1}(a,\eta)}{\partial\mathbf{a}^{2}}\})]. \end{split} \tag{60}$$

In simplified, the above calculation eliminates as described:

*<sup>μ</sup>*2(*α*, *<sup>η</sup>*) = <sup>−</sup>*ξ*(*<sup>n</sup>* <sup>+</sup> *<sup>h</sup>*)*hκ*2*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] *<sup>η</sup><sup>δ</sup>* Γ(*δ* + 1) <sup>+</sup> *<sup>ξ</sup>h*2*κ*4*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] - 2*ξκ*Γ(2*δ* + 1)*η*3*<sup>δ</sup>* (Γ(*<sup>δ</sup>* <sup>+</sup> <sup>1</sup>))2Γ(3*<sup>δ</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> (3 coth2([*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] <sup>−</sup> <sup>1</sup>))*η*2*<sup>δ</sup>* Γ(2*δ* + 1) , *<sup>ν</sup>*2(*α*, *<sup>η</sup>*) = <sup>−</sup>*ξ*(*<sup>n</sup>* <sup>+</sup> *<sup>h</sup>*)*hκ*2*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] coth[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] *<sup>η</sup><sup>δ</sup>* Γ(*δ* + 1) + 1 Γ(*δ* + 1) *<sup>h</sup>*2[2*ξκ*5*cosech*2[*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)]][ *ξκcosech*2(3 coth2([*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] <sup>−</sup> <sup>1</sup>))*η*3*<sup>δ</sup>* Γ(*δ* + 1)Γ(3*δ* + 1) + 2*ξκcosech*<sup>2</sup> coth2([*κ*(*α* + *θ*)])*η*3*<sup>δ</sup>* <sup>Γ</sup>(*<sup>δ</sup>* <sup>+</sup> <sup>1</sup>)Γ(3*<sup>δ</sup>* <sup>+</sup> <sup>1</sup>) <sup>−</sup> <sup>2</sup>*<sup>ξ</sup>* coth(3*cosech*2([*κ*(*<sup>α</sup>* <sup>+</sup> *<sup>θ</sup>*)] <sup>−</sup> <sup>1</sup>))*η*2*<sup>δ</sup>* <sup>Γ</sup>(2*<sup>δ</sup>* <sup>+</sup> <sup>1</sup>) ].

The rest of the iterative terms can be used in the same way. Formerly, the family of q-homotopy analysis transform technique series result of Equation (46) is assumed by:

$$\begin{aligned} \mu(\boldsymbol{a}, \boldsymbol{\eta}) &= \mu\_0(\boldsymbol{a}, \boldsymbol{\eta}) + \sum\_{m=1}^{\infty} \mu\_m(\boldsymbol{a}, \boldsymbol{\eta}) (\frac{1}{n})^m, \\ \nu(\boldsymbol{a}, \boldsymbol{\eta}) &= \nu\_0(\boldsymbol{a}, \boldsymbol{\eta}) + \sum\_{m=1}^{\infty} \nu\_m(\boldsymbol{a}, \boldsymbol{\eta}) (\frac{1}{n})^m. \end{aligned} \tag{61}$$

The exact solution of Equation (46) at *δ* = 1 and taking *ξ* = 0.005, *θ* = 10 and *κ* = 0.1.

$$\begin{aligned} \mu(\mathfrak{a}, \eta) &= = \xi - \kappa \coth[\kappa(\mathfrak{a} + \theta - \xi\eta)], \\ \nu(\mathfrak{a}, \eta) &= -\kappa^2 \cosh^2[\kappa(\mathfrak{a} + \theta - \xi\eta)]. \end{aligned} \tag{62}$$

The solutions *μ*(*α*, *η*) and *ν*(*α*, *η*) are also obtained by using q-HATM and found to be in good agreement with the exact solution of problems. For better understanding the results for both the variables *μ*(*α*, *η*) and *ν*(*α*, *η*) of Example 2 are plotted in Figures 10 and 11 respectively where the higher accuracy is observed.

**Figure 10.** Exact and q-HATM solution of *μ*(*α*, *η*) at *δ* = 1.

**Figure 11.** Exact and q-HATM solution of *ν*(*α*, *η*) at *δ* = 1.

#### **6. Conclusions**

In this paper, we studied the factional view of Whitham–Broer–Kaup equations by using two analytical powerful techniques. With the help of the Laplace and natural transformations, the procedure strengthened and became easy for implementation. A very close contact of the obtained solutions with the exact solution of the problem was observed. It was found that the rate of convergence of the proposed methods was sufficient for solving fractional-order partial differential equations. Therefore, the proposed techniques could be extended to solve other complicated fractional-order problems.

**Author Contributions:** conceptualization, R.S. and H.K.; methodology, R.S.; software, H.K.; validation, D.B., H.K. and R.S.; formal analysis, R.S.; investigation, H.K.; resources, D.B.; data curation, R.S.; writing—original draft preparation, R.S.; writing—review and editing, H.K.; visualization, D.B.; supervision, D.B.

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

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

#### **References**


c 2019 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 (http://creativecommons.org/licenses/by/4.0/).
