New Soliton Solutions of Time-Fractional Korteweg–de Vries Systems

Abstract

Model construction for different physical situations, and developing their solutions, are the major characteristics of the scientific work in physics and engineering. Korteweg–de Vries (KdV) models are very important due to their ability to capture different physical situations such as thin film flows and waves on shallow water surfaces. In this work, a new approach for predicting and analyzing nonlinear time-fractional coupled KdV systems is proposed based on Laplace transform and homotopy perturbation along with Caputo fractional derivatives. This algorithm provides a convergent series solution by applying simple steps through symbolic computations. The efficiency of the proposed algorithm is tested against different nonlinear time-fractional KdV systems, including dispersive long wave and generalized Hirota–Satsuma KdV systems. For validity purposes, the obtained results are compared with the existing solutions from the literature. The convergence of the proposed algorithm over the entire fractional domain is confirmed by finding solutions and errors at various values of fractional parameters. Numerical simulations clearly reassert the supremacy and capability of the proposed technique in terms of accuracy and fewer computations as compared to other available schemes. Analysis reveals that the projected scheme is reliable and hence can be utilized with other kernels in more advanced systems in physics and engineering.

1. Introduction

Kortweg–de Vries (KdV) theory has many applications in physics and engineering. Not surprisingly, many of these applications are in fluid mechanics because many key nonlinear structures were first discovered (chaos, shocks, solitons, bifurcations, etc.), which the underlying theories receive as their most convincing validation [1]. As an integrable nonlinear system, the KdV model has tremendous impact on many parts of theoretical physics, mathematics and the areas in between. Differential equations, algebraic geometry, and Lie group theory have been used on the backhand when dealing with quantum field theory, string theory, and general relatively. KdV systems have also applications in condensed matter, meteorology and plasma physics of protein systems and neuro-physiology [2]. KdV systems are famous for capturing shallow water waves with weak nonlinearity. These equations have many uses in flood analysis and capturing the flows in oceans and vorticity. KdV equations also play a vital role in the study of nonlinear dispersive waves. Initial work on these waves was initiated by Russell and Scott in 1834 [3]; then, Boussinesq and Joseph in the 1870s [4,5,6] explored these waves experimentally. Finally, Korteweg and De Vries re-discovered these waves in 1895 [7]. They first derived them as one-dimensional, long surface gravity and small amplitude waves, which are mostly generated in a shallow water channel. Groesen and Andonowati derived KdV-type of models for water surface waves in [8,9]. Kovalyov and Abadi derived a formula for a class of KdV equations in [10]. The physical behavior of KdV equations is still a hot topic of interest due to its vast applications, and many researchers have attempted to analyze these equations in several phenomena including collision-free hydro magnetic waves, ion-acoustic waves, stratified internal waves, etc. [11,12,13]. One essential KdV model is the Hirota–Satsuma coupled KdV model [14,15] in which the interaction of two long waves has a separate dispersion relationship. It consists of nonlinear PDEs and has various implementations regarding the interaction and evolution of nonlinear waves [16,17,18].

Fractional calculus acts toward the computation of non-integer order derivatives and integrals. Many well-known mathematicians such as Riemann [19], Liouville [20], Euler [21], Laurent [22], Laplace [23], etc. have contributed a lot in this field. Fractional calculus plays an important role in different areas of science and engineering such as diffusion procedures, decentralized wireless networks, biology, physics, viscoelasticity, bio-engineering, micro-grids, geochemistry, etc. Most of the mentioned physical models ultimately convert to fractional order systems of differential equations. Solving such systems for prediction and analysis purposes is a very important area of research. Many researchers have introduced several methods for the solution of fractional order systems. The variation iteration method (VIM) [24], Adomian decomposition method (ADM) [25], auxiliary equation method [26], homotopy analysis method (HAM) [27], and homotopy perturbation method (HPM) [28] are some of these methods. Wang attempted time and spatial fractional KdV–Burgers equation through HPM in [29]. Kurulay et al. a studied fractional modified KdV using the differential transform method in [30]. El-Ajou et al. investigated a nonlinear space-time-fractional KdV–Burgers equation using a residual power series algorithm in [31]. Elbadri et al. used natural decomposition for the solution of time-fractional KdV systems in [32]. Mohammed et al. used the Riccati Equation Method for the solution of a stochastic fractional Kuramoto–Sivashinsky model in [33].

HPM is the combination of homotopy with a perturbation process which provides a better way to find series solutions of differential equations. Several modifications of classical HPM are introduced in the literature with the aim to improve accuracy [34,35,36]. One of the modifications is the Laplace Homotopy Perturbation Method (LHPM) [37], in which HPM is hybridized with Laplace transformation. In this continuation, LHPM is extended to nonlinear coupled time-fractional KdV systems for better results in terms of accuracy with less computations. In the rest of the manuscript, Section 2 and Section 3 present preliminaries and the basic theory of LHPM for fractional-order systems, respectively. Section 4 provides an application of LHPM to nonlinear time-fractional KDV systems. The discussion and conclusion are given in Section 5 and Section 6, respectively.

2. Preliminaries

Definition 1.

The Laplace transform $\mathcal{L}\left[\mathcal{R}\right(\mathfrak{u},\mathfrak{t}\left)\right]$ of Riemann–Liouville time-fractional integral ${\mathcal{I}}_t^{\beta }$ is defined as follows [38]:

$$\mathcal{L}\left[{\mathcal{I}}_t^{\beta }\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]=s^{-\beta }\mathcal{L}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]$$

(1)

Definition 2.

The Laplace transform $\mathcal{L}\left[\mathcal{R}\right(\mathfrak{u},\mathfrak{t}\left)\right]$ of Caputo’s time-fractional derivative ${\mathcal{D}}_t^{\beta }$ is defined as follows [38]:

$$\mathcal{L}\left[{\mathcal{D}}_t^{\beta }\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]=s^{\beta }\mathcal{L}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]-\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{R}}^{\left(q\right)}(\mathfrak{u},0),q-1<\beta \le q.$$

(2)

3. Basic Idea of LHPM for Time-Fractional Systems

In order to understand the basic concept of the proposed algorithm, let us consider the general time-fractional system of partial differential equations as:

$$\begin{array}{cc}\hfill {\mathcal{D}}_t^{\beta }\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{L}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{h}(\mathfrak{u},\mathfrak{t})& =0,\hfill \\ \hfill {\mathcal{D}}_t^{\beta }\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{L}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{g}(\mathfrak{u},\mathfrak{t})& =0,\hfill \\ \hfill {\mathcal{D}}_t^{\beta }\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{L}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{f}(\mathfrak{u},\mathfrak{t})& =0,\mathfrak{u}ϵ\Omega ,\mathfrak{t}>0,q-1<\beta \le q,\hfill \end{array}$$

(3)

with initial conditions

$$\begin{array}{c}\hfill \mathcal{R}(\mathfrak{u},0)=\mathcal{I}1,\\ \hfill \mathcal{S}(\mathfrak{u},0)=\mathcal{I}2,\\ \hfill \mathcal{T}(\mathfrak{u},0)=\mathcal{I}3.\end{array}$$

(4)

where $\mathcal{R}(\mathfrak{u},\mathfrak{t})$, $\mathcal{S}(\mathfrak{u},\mathfrak{t})$ and $\mathcal{T}(\mathfrak{u},\mathfrak{t})$ are unknown functions and ${\mathcal{D}}_t^{\beta }$ is the fractional-order derivative, $\mathfrak{h}(\mathfrak{u},\mathfrak{t})$, $\mathfrak{g}(\mathfrak{u},\mathfrak{t})$ and $\mathfrak{f}(\mathfrak{u},\mathfrak{t})$ are known functions, $\mathcal{I}i,i=1,2,3$ are initial conditions, and $\Omega$ is the spatial domain. In addition, $\mathcal{L}$ and $\mathcal{N}$ represent the linear and nonlinear operators, respectively.

The first step of the algorithm is to apply Laplace transformation to the system (3) as

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathcal{D}}_t^{\beta }\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]\right\}+\mathcal{L}\{\mathcal{L}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{h}(\mathfrak{u},\mathfrak{t})\}=0,\\ \hfill \mathcal{L}\left\{{\mathcal{D}}_t^{\beta }\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]\right\}+\mathcal{L}\{\mathcal{L}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{g}(\mathfrak{u},\mathfrak{t})\}=0,\\ \hfill \mathcal{L}\left\{{\mathcal{D}}_t^{\beta }\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]\right\}+\mathcal{L}\{\mathcal{L}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{f}(\mathfrak{u},\mathfrak{t})\}=0,\end{array}$$

(5)

By using definition (2), we will transform the following system in the s-domain

$$\begin{array}{c}\hfill \mathcal{L}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{R}}^{\left(q\right)}(\mathfrak{u},0)+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{h}(\mathfrak{u},\mathfrak{t})\}=0,\\ \hfill \mathcal{L}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{S}}^{\left(q\right)}(\mathfrak{u},0)+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{g}(\mathfrak{u},\mathfrak{t})\}=0,\\ \hfill \mathcal{L}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{T}}^{\left(q\right)}(\mathfrak{u},0)+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{f}(\mathfrak{u},\mathfrak{t})\}=0,\end{array}$$

(6)

In next step of the algorithm, let us construct a set of homotopies as

$$\begin{array}{c}\hfill {\mathcal{H}}_1=(1-p)(\mathcal{L}\left\{\mathcal{R}(\mathfrak{u},\mathfrak{t};p)\right\}-{\mathcal{R}}_0(\mathfrak{u},\mathfrak{t}))+p(\mathcal{L}\left\{\mathcal{R}(\mathfrak{u},\mathfrak{t};p)\right\}-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{R}}^{\left(q\right)}(\mathfrak{u},0)+\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{R}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{h}(\mathfrak{u},\mathfrak{t})\}),\\ \hfill {\mathcal{H}}_2=(1-p)(\mathcal{L}\left\{\mathcal{S}(\mathfrak{u},\mathfrak{t};p)\right\}-{\mathcal{S}}_0(\mathfrak{u},\mathfrak{t}))+p(\mathcal{L}\left\{\mathcal{S}(\mathfrak{u},\mathfrak{t};p)\right\}-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{S}}^{\left(q\right)}(\mathfrak{u},0)+\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{S}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{g}(\mathfrak{u},\mathfrak{t})\}),\\ \hfill {\mathcal{H}}_3=(1-p)(\mathcal{L}\left\{\mathcal{T}(\mathfrak{u},\mathfrak{t};p)\right\}-{\mathcal{T}}_0(\mathfrak{u},\mathfrak{t}))+p(\mathcal{L}\left\{\mathcal{T}(\mathfrak{u},\mathfrak{t};p)\right\}-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{T}}^{\left(q\right)}(\mathfrak{u},0)+\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[\mathcal{T}(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{f}(\mathfrak{u},\mathfrak{t})\}).\end{array}$$

(7)

where ${\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})$, ${\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})$ and ${\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})$ are initial guesses that satisfy the given conditions. For $p=0$ in (7), we have

$$\begin{array}{c}\hfill \mathcal{L}\left\{\mathcal{R}(\mathfrak{u},\mathfrak{t};0)\right\}={\mathcal{R}}_0(\mathfrak{u},\mathfrak{t}),\\ \hfill \mathcal{L}\left\{\mathcal{S}(\mathfrak{u},\mathfrak{t};0)\right\}={\mathcal{S}}_0(\mathfrak{u},\mathfrak{t}),\\ \hfill \mathcal{L}\left\{\mathcal{T}(\mathfrak{u},\mathfrak{t};0)\right\}={\mathcal{T}}_0(\mathfrak{u},\mathfrak{t}).\end{array}$$

(8)

and $p=1$ in (7) gives

$$\begin{array}{c}\hfill \mathcal{L}\left\{\mathcal{R}\right(\mathfrak{u},\mathfrak{t};1\left)\right\}=\mathcal{R}(\mathfrak{u},\mathfrak{t}),\\ \hfill \mathcal{L}\left\{\mathcal{S}\right(\mathfrak{u},\mathfrak{t};1\left)\right\}=\mathcal{S}(\mathfrak{u},\mathfrak{t}),\\ \hfill \mathcal{L}\left\{\mathcal{T}\right(\mathfrak{u},\mathfrak{t};1\left)\right\}=\mathcal{T}(\mathfrak{u},\mathfrak{t}).\end{array}$$

(9)

By expanding $\mathcal{R}(\mathfrak{u},\mathfrak{t};p)$, $\mathcal{S}(\mathfrak{u},\mathfrak{t};p)$ and $\mathcal{T}(\mathfrak{u},\mathfrak{t};p)$ as Taylor series with respect to p, we have

$$\begin{array}{c}\hfill \mathcal{R}(\mathfrak{u},\mathfrak{t};p)=\sum _{m=1}^{\infty }{p}^m{\mathcal{R}}_m,\\ \hfill \mathcal{S}(\mathfrak{u},\mathfrak{t};p)=\sum _{m=1}^{\infty }{p}^m{\mathcal{S}}_m,\\ \hfill \mathcal{T}(\mathfrak{u},\mathfrak{t};p)=\sum _{m=1}^{\infty }{p}^m{\mathcal{T}}_m.\end{array}$$

(10)

Next, substituting (10) in (7) and comparing like powers of p will give the following deformations.

The first-order deformation is

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathcal{R}}_1(\mathfrak{u},\mathfrak{t})\right\}+{\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{R}}^{\left(q\right)}(\mathfrak{u},0)+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})\right]+\\ \hfill \mathcal{N}\left[{\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{h}(\mathfrak{u},\mathfrak{t})\}=0,\\ \hfill \mathcal{L}\left\{{\mathcal{S}}_1(\mathfrak{u},\mathfrak{t})\right\}+{\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{S}}^{\left(q\right)}(\mathfrak{u},0)+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})\right]+\\ \hfill \mathcal{N}\left[{\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{g}(\mathfrak{u},\mathfrak{t})\}=0,\\ \hfill \mathcal{L}\left\{{\mathcal{T}}_1(\mathfrak{u},\mathfrak{t})\right\}+{\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})-\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{T}}^{\left(q\right)}(\mathfrak{u},0)+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})\right]+\\ \hfill \mathcal{N}\left[{\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{f}(\mathfrak{u},\mathfrak{t})\}=0.\end{array}$$

(11)

An application of inverse Laplace transform leads to

$$\begin{array}{c}\hfill {\mathcal{R}}_1(\mathfrak{u},\mathfrak{t})={\mathcal{L}}^{-1}\left\{-{\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})+\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{R}}^{\left(q\right)}(\mathfrak{u},0)\right\}-\\ \hfill {\mathcal{L}}^{-1}\left\{\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{h}(\mathfrak{u},\mathfrak{t})\right\},\\ \hfill {\mathcal{S}}_1(\mathfrak{u},\mathfrak{t})={\mathcal{L}}^{-1}\left\{-{\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})+\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{S}}^{\left(q\right)}(\mathfrak{u},0)\right\}-\\ \hfill {\mathcal{L}}^{-1}\left\{\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{g}(\mathfrak{u},\mathfrak{t})\right\},\\ \hfill {\mathcal{T}}_1(\mathfrak{u},\mathfrak{t})={\mathcal{L}}^{-1}\left\{-{\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})+\left(\frac{1}{s^{\beta }}\right)\sum _{p=0}^{q-1}{s}^{\beta -p-1}{\mathcal{T}}^{\left(q\right)}(\mathfrak{u},0)\right\}-\\ \hfill {\mathcal{L}}^{-1}\left\{\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})\right]-\mathfrak{f}(\mathfrak{u},\mathfrak{t})\right\}.\end{array}$$

(12)

Similarly, a kth-order deformation is

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathcal{R}}_k(\mathfrak{u},\mathfrak{t})\right\}+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{R}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{R}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]\}=0,\\ \hfill \mathcal{L}\left\{{\mathcal{S}}_k(\mathfrak{u},\mathfrak{t})\right\}+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{S}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{S}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]\}=0,\\ \hfill \mathcal{L}\left\{{\mathcal{T}}_k(\mathfrak{u},\mathfrak{t})\right\}+\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{T}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{T}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]\}=0.\end{array}$$

(13)

Operating inverse Laplace transform gives the following

$$\begin{array}{cc}\hfill {\mathcal{R}}_k(\mathfrak{u},\mathfrak{t})& ={\mathcal{L}}^{-1}\left\{-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{R}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{R}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]\}\right\},\hfill \\ \hfill {\mathcal{S}}_k(\mathfrak{u},\mathfrak{t})& ={\mathcal{L}}^{-1}\left\{-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{S}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{S}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]\}\right\}.\hfill \\ \hfill {\mathcal{T}}_k(\mathfrak{u},\mathfrak{t})& ={\mathcal{L}}^{-1}\left\{-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\mathcal{L}\left[{\mathcal{T}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]+\mathcal{N}\left[{\mathcal{T}}_{k-1}(\mathfrak{u},\mathfrak{t})\right]\}\right\}.\hfill \end{array}$$

(14)

The approximate series solution of the fractional system is

$$\begin{array}{cc}\hfill \tilde{\mathcal{R}}& ={\mathcal{R}}_0(\mathfrak{u},\mathfrak{t})+{\mathcal{R}}_1(\mathfrak{u},\mathfrak{t})+{\mathcal{R}}_2(\mathfrak{u},\mathfrak{t})+{\mathcal{R}}_3(\mathfrak{u},\mathfrak{t})+\cdots \hfill \\ \hfill \tilde{\mathcal{S}}& ={\mathcal{S}}_0(\mathfrak{u},\mathfrak{t})+{\mathcal{S}}_1(\mathfrak{u},\mathfrak{t})+{\mathcal{S}}_2(\mathfrak{u},\mathfrak{t})+{\mathcal{S}}_3(\mathfrak{u},\mathfrak{t})+\cdots \hfill \\ \hfill \tilde{\mathcal{T}}& ={\mathcal{T}}_0(\mathfrak{u},\mathfrak{t})+{\mathcal{T}}_1(\mathfrak{u},\mathfrak{t})+{\mathcal{T}}_2(\mathfrak{u},\mathfrak{t})+{\mathcal{T}}_3(\mathfrak{u},\mathfrak{t})+\cdots \hfill \end{array}$$

(15)

Residual errors can be found by substituting approximate solutions back in the given system as

$$\begin{array}{cc}\hfill \mathfrak{Res}1& ={\mathcal{D}}_t^{\beta }\left[\tilde{\mathcal{R}}\right]+\mathcal{L}\left[\tilde{\mathcal{R}}\right]+\mathcal{N}\left[\tilde{\mathcal{R}}\right]-\mathfrak{h}(\mathfrak{u},\mathfrak{t}),\hfill \\ \hfill \mathfrak{Res}2& ={\mathcal{D}}_t^{\beta }\left[\tilde{\mathcal{S}}\right]+\mathcal{L}\left[\tilde{\mathcal{S}}\right]+\mathcal{N}\left[\tilde{\mathcal{S}}\right]-\mathfrak{g}(\mathfrak{u},\mathfrak{t}),\hfill \\ \hfill \mathfrak{Res}3& ={\mathcal{D}}_t^{\beta }\left[\tilde{\mathcal{T}}\right]+\mathcal{L}\left[\tilde{\mathcal{T}}\right]+\mathcal{N}\left[\tilde{\mathcal{T}}\right]-\mathfrak{f}(\mathfrak{u},\mathfrak{t}).\hfill \end{array}$$

(16)

Similarly, we can extend the proposed algorithm to a system containing more than three equations.

4. Application LHPM to Fractional KdV Systems

Example 1.

Consider the time-fractional, nonlinear coupled KdV system

$$\begin{array}{cc}\hfill \frac{{\partial }^{\beta }\mathbf{U}}{\partial {\mathfrak{t}}^{\beta }}& =-\mathfrak{a}\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-6\mathfrak{a}\mathbf{U}\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}+6\mathbf{V}\frac{\partial \mathbf{V}}{\partial \mathfrak{x}},\hfill \\ \hfill \frac{{\partial }^{\beta }\mathbf{V}}{\partial {\mathfrak{t}}^{\beta }}& =-\mathfrak{a}\frac{{\partial }^3\mathbf{V}}{\partial {\mathfrak{x}}^3}-3\mathfrak{a}\mathbf{U}\frac{\partial \mathbf{V}}{\partial \mathfrak{x}},0<\beta \le 1,\hfill \end{array}$$

(17)

subject to initial conditions

$$\begin{array}{cc}\hfill \mathbf{U}(\mathfrak{x},0)& ={\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right),\hfill \\ \hfill \mathbf{V}(\mathfrak{x},0)& =\sqrt{\frac{\mathfrak{a}}{2}}{\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right),\hfill \end{array}$$

(18)

where α and η are arbitrary constants and $\mathfrak{a}$ is wave velocity. For β = 1, the exact solution of (17) is

$$\begin{array}{cc}\hfill \mathbf{U}(\mathfrak{x},\mathfrak{t})& ={\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}-\frac{\mathfrak{a}{\eta }^3\mathfrak{t}}{2}\right),\hfill \\ \hfill \mathbf{V}(\mathfrak{x},\mathfrak{t})& =\sqrt{\frac{\mathfrak{a}}{2}}{\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}-\frac{\mathfrak{a}{\eta }^3\mathfrak{t}}{2}\right).\hfill \end{array}$$

(19)

Solution: Firstly, apply Laplace transform to a system (17) and use definition (2) to obtain the following

$$\begin{array}{c}\hfill s^{\beta }\mathcal{L}\left[\mathbf{U}(\mathfrak{x},\mathfrak{t})\right]-s^{\beta -1}{\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)=\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-6\mathfrak{a}\mathbf{U}\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}+6\mathbf{V}\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}\right\},\\ \hfill s^{\beta }\mathcal{L}\left[\mathbf{V}(\mathfrak{x},\mathfrak{t})\right]-s^{\beta -1}\sqrt{\frac{\mathfrak{a}}{2}}{\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)=\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3\mathbf{V}}{\partial {\mathfrak{x}}^3}-3\mathfrak{a}\mathbf{U}\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}\right\},\end{array}$$

(20)

In the next step, construct a set of homotopies for the system as

$$\begin{array}{c}\hfill {\mathcal{H}}_1=(1-p)(\mathcal{L}\left\{\mathbf{U}(\mathfrak{x},\mathfrak{t})\right\}-{\mathbf{U}}_0(\mathfrak{x},\mathfrak{t}))+p(\mathcal{L}\left\{\mathbf{U}(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s}\right){\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)-\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-6\mathfrak{a}\mathbf{U}\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}+6\mathbf{V}\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}\right\}).\end{array}$$

(21)

$$\begin{array}{c}\hfill {\mathcal{H}}_2=(1-p)(\mathcal{L}\left\{\mathbf{V}(\mathfrak{x},\mathfrak{t})\right\}-{\mathbf{V}}_0(\mathfrak{x},\mathfrak{t}))+p(\mathcal{L}\left\{\mathbf{V}(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s}\right)\sqrt{\frac{\mathfrak{a}}{2}}{\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)-\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3\mathbf{V}}{\partial {\mathfrak{x}}^3}-3\mathfrak{a}\mathbf{U}\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}\right\}),\end{array}$$

(22)

where

$$\begin{array}{cc}\hfill {\mathbf{U}}_0(\mathfrak{x},\mathfrak{t})& ={\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right),\hfill \\ \hfill {\mathbf{V}}_0(\mathfrak{x},\mathfrak{t})& =\sqrt{\frac{\mathfrak{a}}{2}}{\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right),\hfill \end{array}$$

(23)

are the initial guesses which satisfies (18). Using (10) in (21) and (22) and then comparing coefficients for same power of p will lead to the following various order problems.

First-order deformation is

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_1(\mathfrak{x},\mathfrak{t})\right\}+{\mathbf{U}}_0(\mathfrak{x},\mathfrak{t})-\left(\frac{1}{s}\right){\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{-\mathfrak{a}\frac{{\partial }^3{\mathbf{U}}_0}{\partial {\mathfrak{x}}^3}-\\ \hfill 6\mathfrak{a}{\mathbf{U}}_0\frac{\partial {\mathbf{U}}_0}{\partial \mathfrak{x}}+6{\mathbf{V}}_0\frac{\partial {\mathbf{V}}_0}{\partial \mathfrak{x}}\}=0,\end{array}$$

(24)

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{V}}_1(\mathfrak{x},\mathfrak{t})\right\}+{\mathbf{V}}_0(\mathfrak{x},\mathfrak{t})-\left(\frac{1}{s}\right)\sqrt{\frac{\mathfrak{a}}{2}}{\eta }^{2}sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{-\mathfrak{a}\frac{{\partial }^3{\mathbf{V}}_0}{\partial {\mathfrak{x}}^3}-\\ \hfill 3\mathfrak{a}{\mathbf{U}}_0\frac{\partial {\mathbf{V}}_0}{\partial \mathfrak{x}}\}=0,\end{array}$$

(25)

with conditions

$$\begin{array}{cc}\hfill {\mathbf{U}}_1(\mathfrak{x},0)& =0,\hfill \\ \hfill {\mathbf{V}}_1(\mathfrak{x},0)& =0,\hfill \end{array}$$

(26)

Application of inverse Laplace transformation will give the following

$$\begin{array}{cc}\hfill {\mathbf{U}}_1(\mathfrak{x},\mathfrak{t})& =\frac{{\mathfrak{t}}^{\beta }\mathfrak{a}{\eta }^5\left(tanh\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)sec{h}^4\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)+tan{h}^3\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)\right)}{\Gamma (\beta +1)},\hfill \\ \hfill {\mathbf{V}}_1(\mathfrak{x},\mathfrak{t})& =\frac{{\mathfrak{t}}^{\beta }{\mathfrak{a}}^{3/2}{\eta }^5\left(tanh\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)sec{h}^4\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)+tan{h}^3\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)sec{h}^2\left(\frac{\alpha }{2}+\frac{\eta \mathfrak{x}}{2}\right)\right)}{\sqrt{2}\Gamma (\beta +1)},\hfill \end{array}$$

(27)

Second-order deformation is

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_2(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3{\mathbf{U}}_1}{\partial {\mathfrak{x}}^3}-6\mathfrak{a}{\mathbf{U}}_1\frac{\partial {\mathbf{U}}_1}{\partial \mathfrak{x}}+6{\mathbf{V}}_1\frac{\partial {\mathbf{V}}_1}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_2(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3{\mathbf{V}}_1}{\partial {\mathfrak{x}}^3}-3\mathfrak{a}{\mathbf{U}}_1\frac{\partial {\mathbf{V}}_1}{\partial \mathfrak{x}}\right\}=0,\\ \hfill {\mathbf{U}}_2(\mathfrak{x},0)=0,{\mathbf{V}}_2(\mathfrak{x},0)=0.\end{array}$$

(28)

Application of inverse Laplace transformation will give the following solution component

$$\begin{array}{cc}\hfill {\mathbf{U}}_2(\mathfrak{x},\mathfrak{t})& =\frac{{\mathfrak{t}}^{2\beta }{\mathfrak{a}}^2{\eta }^8(cosh(\alpha +\eta \mathfrak{x})-2)sec{h}^4\left(\frac{\alpha +\eta \mathfrak{x}}{2}\right)}{2\Gamma (2\beta +1)},\hfill \\ \hfill {\mathbf{V}}_2(\mathfrak{x},\mathfrak{t})& =\frac{{\mathfrak{t}}^{2\beta }{\mathfrak{a}}^{5/2}{\eta }^8(cosh(\alpha +\eta \mathfrak{x})-2)sec{h}^4\left(\frac{\alpha +\eta \mathfrak{x}}{2}\right)}{2\sqrt{2}\Gamma (2\beta +1)},\hfill \end{array}$$

(29)

Third-order deformation is

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_3(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3{\mathbf{U}}_2}{\partial {\mathfrak{x}}^3}-6\mathfrak{a}{\mathbf{U}}_2\frac{\partial {\mathbf{U}}_2}{\partial \mathfrak{x}}+6{\mathbf{V}}_2\frac{\partial {\mathbf{V}}_2}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_3(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\mathfrak{a}\frac{{\partial }^3{\mathbf{V}}_2}{\partial {\mathfrak{x}}^3}-3\mathfrak{a}{\mathbf{U}}_2\frac{\partial {\mathbf{V}}_2}{\partial \mathfrak{x}}\right\}=0,\\ \hfill {\mathbf{U}}_3(\mathfrak{x},0)=0,{\mathbf{V}}_3(\mathfrak{x},0)=0.\end{array}$$

(30)

Application of inverse Laplace transformation will give the following

$$\begin{array}{c}\hfill {\mathbf{U}}_3(\mathfrak{x},\mathfrak{t})={\mathfrak{a}}^3{\eta }^{11}{\mathfrak{t}}^{3\beta }tanh\left(\frac{\alpha +\eta \mathfrak{x}}{2}\right)sec{h}^6\left(\frac{\alpha +\eta \mathfrak{x}}{2}\right)(\Gamma {(\beta +1)}^2(-32cosh(\alpha +\eta \mathfrak{x})+\\ \hfill {cosh\left(2(\alpha +\eta \mathfrak{x})\right)+39)+12\Gamma (2\beta +1)(cosh(\alpha +\eta \mathfrak{x})-2))/8\Gamma (\beta +1)}^2\Gamma (3\beta +1),\\ \hfill {\mathbf{V}}_3(\mathfrak{x},\mathfrak{t})={\mathfrak{a}}^{7/2}{\eta }^{11}{\mathfrak{t}}^{3\beta }tanh\left(\frac{\alpha +\eta \mathfrak{x}}{2}\right)sec{h}^6\left(\frac{\alpha +\eta \mathfrak{x}}{2}\right)(\Gamma {(\beta +1)}^2(-32cosh(\alpha +\eta \mathfrak{x})+\\ \hfill cosh\left(2(\alpha +\eta \mathfrak{x})\right)+39)+12\Gamma (2\beta +1)(cosh(\alpha +\eta \mathfrak{x})-2))/8\sqrt{2}\Gamma {(\beta +1)}^2\Gamma (3\beta +1),\end{array}$$

(31)

Proceeding similarly, higher-order problems and solutions can be obtained. For the current manuscript, we approximate all the solutions up to the fifth order.

$$\begin{array}{c}\hfill \tilde{\mathbf{U}}=\sum _{m=0}^5{\mathbf{U}}_m(\mathfrak{x},\mathfrak{t}),\\ \hfill \tilde{\mathbf{V}}=\sum _{m=0}^5{\mathbf{V}}_m(\mathfrak{x},\mathfrak{t}).\end{array}$$

(32)

Residual errors at the fifth order can be obtained by substituting the approximate solution in (17)

$$\begin{array}{cc}\hfill \mathfrak{Res}\mathbf{U}& =-\mathfrak{a}\frac{{\partial }^3\tilde{\mathbf{U}}}{\partial {\mathfrak{x}}^3}-6\mathfrak{a}\tilde{\mathbf{U}}\frac{\partial \tilde{\mathbf{U}}}{\partial \mathfrak{x}}+6\tilde{\mathbf{V}}\frac{\partial \tilde{\mathbf{V}}}{\partial \mathfrak{x}},\hfill \\ \hfill \mathfrak{Res}\mathbf{V}& =-\mathfrak{a}\frac{{\partial }^3\tilde{\mathbf{V}}}{\partial {\mathfrak{x}}^3}-3\mathfrak{a}\tilde{\mathbf{U}}\frac{\partial \tilde{\mathbf{V}}}{\partial \mathfrak{x}}.\hfill \end{array}$$

(33)

Example 2.

Consider a time-fractional, nonlinear dispersive long-wave system

$$\begin{array}{cc}\hfill \frac{{\partial }^{\beta }\mathbf{U}}{\partial {\mathfrak{t}}^{\beta }}& =-\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}-\frac{1}{2}\frac{\partial {\mathbf{U}}^2}{\partial \mathfrak{x}},\hfill \\ \hfill \frac{{\partial }^{\beta }\mathbf{V}}{\partial {\mathfrak{t}}^{\beta }}& =-\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}-\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-\frac{\partial \mathbf{UV}}{\partial \mathfrak{x}},0<\beta \le 1,\hfill \end{array}$$

(34)

subject to initial conditions

$$\begin{array}{cc}\hfill \mathbf{U}(\mathfrak{x},0)& =\mathfrak{a}tanh\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)+\mathfrak{a},\hfill \\ \hfill \mathbf{V}(\mathfrak{x},0)& =-1+\frac{1}{2}{\mathfrak{a}}^{2}sec{h}^2\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right),\hfill \end{array}$$

(35)

where $\mathfrak{a}$ and η are arbitrary constants. For β = 1, the exact solution of (34) is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{U}(\mathfrak{x},\mathfrak{t})& =\mathfrak{a}tanh\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}-\frac{{\mathfrak{a}}^2\mathfrak{t}}{2}\right)+\mathfrak{a},\hfill \\ \hfill \mathbf{V}(\mathfrak{x},\mathfrak{t})& =-1+\frac{1}{2}{\mathfrak{a}}^{2}sec{h}^2\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}-\frac{{\mathfrak{a}}^2\mathfrak{t}}{2}\right).\hfill \end{array}\end{array}$$

(36)

Solution: Application of Laplace transformation on both sides of system (34) as

$$\begin{array}{cc}\hfill \mathcal{L}\left\{\frac{{\partial }^{\beta }\mathbf{U}}{\partial {\mathfrak{t}}^{\beta }}\right\}& =\mathcal{L}\left\{-\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}-\frac{1}{2}\frac{\partial {\mathbf{U}}^2}{\partial \mathfrak{x}}\right\},\hfill \\ \hfill \mathcal{L}\left\{\frac{{\partial }^{\beta }\mathbf{V}}{\partial {\mathfrak{t}}^{\beta }}\right\}& =\mathcal{L}\left\{-\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}-\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-\frac{\partial \mathbf{UV}}{\partial \mathfrak{x}}\right\},\hfill \end{array}$$

(37)

After using (2), (37) is reduced to the following

$$\begin{array}{c}\hfill s^{\beta }\mathcal{L}\left[\mathbf{U}(\mathfrak{x},\mathfrak{t})\right]-s^{\beta -1}\mathfrak{a}tanh\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)+\mathfrak{a}=\mathcal{L}\left\{-\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}-\frac{1}{2}\frac{\partial {\mathbf{U}}^2}{\partial \mathfrak{x}}\right\},\\ \hfill s^{\beta }\mathcal{L}\left[\mathbf{V}(\mathfrak{x},\mathfrak{t})\right]-s^{\beta -1}-1+\frac{1}{2}{\mathfrak{a}}^{2}sec{h}^2\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)=\mathcal{L}\left\{-\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}-\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-\frac{\partial \mathbf{UV}}{\partial \mathfrak{x}}\right\},\end{array}$$

(38)

Next, we need to construct a set of homotopies for the fractional system as

$$\begin{array}{c}\hfill {\mathcal{H}}_1=(1-p)(\mathcal{L}\left\{\mathbf{U}(\mathfrak{x},\mathfrak{t})\right\}-{\mathbf{U}}_0(\mathfrak{x},\mathfrak{t}))+p(\mathcal{L}\left\{\mathbf{U}(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s}\right)\left(\mathfrak{a}tanh\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)+\mathfrak{a}\right)-\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{\partial \mathbf{V}}{\partial \mathfrak{x}}-\frac{1}{2}\frac{\partial {\mathbf{U}}^2}{\partial \mathfrak{x}}\right\}),\end{array}$$

(39)

$$\begin{array}{c}\hfill {\mathcal{H}}_2=(1-p)(\mathcal{L}\left\{\mathbf{V}(\mathfrak{x},\mathfrak{t})\right\}-{\mathbf{V}}_0(\mathfrak{x},\mathfrak{t}))+p(\mathcal{L}\left\{\mathbf{V}(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s}\right)\left(-1+\frac{1}{2}{\mathfrak{a}}^{2}sec{h}^2\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)\right)-\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}-\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-\frac{\partial \mathbf{UV}}{\partial \mathfrak{x}}\right\}),\end{array}$$

(40)

where

$$\begin{array}{cc}\hfill {\mathbf{U}}_0(\mathfrak{x},\mathfrak{t})& =\mathfrak{a}tanh\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)+\mathfrak{a},\hfill \\ \hfill {\mathbf{V}}_0(\mathfrak{x},\mathfrak{t})& =-1+\frac{1}{2}{\mathfrak{a}}^{2}sec{h}^2\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right),\hfill \end{array}$$

(41)

are the initial guesses for the current problem.

In the next step of the algorithm, by using (10) in (40) and (39), then comparing coefficients for same power of p, we obtain the following different order problems.

First-order problem:

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_1(\mathfrak{x},\mathfrak{t})\right\}+{\mathbf{U}}_0(\mathfrak{x},\mathfrak{t})-\left(\frac{1}{s}\right)\left(\mathfrak{a}tanh\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)+\mathfrak{a}\right)-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{-\frac{\partial {\mathbf{V}}_0}{\partial \mathfrak{x}}-\\ \hfill \frac{1}{2}\frac{\partial {\mathbf{U}}_0^2}{\partial \mathfrak{x}}\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_1(\mathfrak{x},\mathfrak{t})\right\}+{\mathbf{V}}_0(\mathfrak{x},\mathfrak{t})-\left(\frac{1}{s}\right)\left(-1+\frac{1}{2}{\mathfrak{a}}^{2}sec{h}^2\left(\frac{\eta }{2}+\frac{\mathfrak{ax}}{2}\right)\right)-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{-\frac{\partial {\mathbf{U}}_0}{\partial \mathfrak{x}}\\ \hfill -\frac{{\partial }^3{\mathbf{U}}_0}{\partial {\mathfrak{x}}^3}-\frac{\partial {\mathbf{U}}_0{\mathbf{V}}_0}{\partial \mathfrak{x}}\}=0,\\ \hfill {\mathbf{U}}_1(\mathfrak{x},0)=0,{\mathbf{V}}_1(\mathfrak{x},0)=0.\end{array}$$

(42)

Second-order problem:

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_2(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{\partial {\mathbf{V}}_1}{\partial \mathfrak{x}}-\frac{1}{2}\frac{\partial {\mathbf{U}}_1^2}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_2(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{\partial {\mathbf{U}}_1}{\partial \mathfrak{x}}-\frac{{\partial }^3{\mathbf{U}}_1}{\partial {\mathfrak{x}}^3}-\frac{\partial {\mathbf{U}}_1{\mathbf{V}}_1}{\partial \mathfrak{x}}\right\}=0,\\ \hfill {\mathbf{U}}_2(\mathfrak{x},0)=0,{\mathbf{V}}_2(\mathfrak{x},0)=0.\end{array}$$

(43)

Third-order problem:

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_3(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{\partial {\mathbf{V}}_2}{\partial \mathfrak{x}}-\frac{1}{2}\frac{\partial {\mathbf{U}}_2^2}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_3(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{\partial {\mathbf{U}}_2}{\partial \mathfrak{x}}-\frac{{\partial }^3{\mathbf{U}}_2}{\partial {\mathfrak{x}}^3}-\frac{\partial {\mathbf{U}}_2{\mathbf{V}}_2}{\partial \mathfrak{x}}\right\}=0,\\ \hfill {\mathbf{U}}_3(\mathfrak{x},0)=0,{\mathbf{V}}_3(\mathfrak{x},0)=0.\end{array}$$

(44)

Continuing this way, we can find higher-order problems. The approximate series solution of fifth order is obtained in this example, and numerical results are shown in Tables 3 and 4 and Figures 4–6.

Example 3.

Consider a time-fractional generalized Hirota–Satsuma coupled KdV system

$$\begin{array}{cc}\hfill \frac{{\partial }^{\beta }\mathbf{U}}{\partial {\mathfrak{t}}^{\beta }}& =\frac{1}{2}\frac{{\partial }^3\mathbf{U}}{\partial {\mathfrak{x}}^3}-3\mathbf{U}\frac{\partial \mathbf{U}}{\partial \mathfrak{x}}+3\frac{\partial \left(\mathbf{V}\mathbf{W}\right)}{\partial \mathfrak{x}},\hfill \\ \hfill \frac{{\partial }^{\beta }\mathbf{V}}{\partial {\mathfrak{t}}^{\beta }}& =-\frac{{\partial }^3\mathbf{V}}{\partial {\mathfrak{x}}^3}+3\mathbf{U}\frac{\partial \mathbf{V}}{\partial \mathfrak{x}},\hfill \\ \hfill \frac{{\partial }^{\beta }\mathbf{W}}{\partial {\mathfrak{t}}^{\beta }}& =-\frac{{\partial }^3\mathbf{W}}{\partial {\mathfrak{x}}^3}+3\mathbf{U}\frac{\partial \mathbf{W}}{\partial \mathfrak{x}},0<\beta \le 1,\hfill \end{array}$$

(45)

subject to initial conditions

$$\begin{array}{cc}\hfill \mathbf{U}(\mathfrak{x},0)& =\frac{\mathfrak{p}-8{\alpha }^2}{3}+4{\alpha }^{2}tan{h}^2\left(\alpha \mathfrak{x}\right),\hfill \\ \hfill \mathbf{V}(\mathfrak{x},0)& =-\frac{4{\alpha }^2(3{\alpha }^2{\mathfrak{c}}_0-2\mathfrak{p}{\mathfrak{c}}_2+4{\alpha }^2{\mathfrak{c}}_2)}{3{\mathfrak{c}}_2^2}+{\frac{4{\alpha }^4}{\mathfrak{c}}}_{2}tan{h}^2\left(\alpha \mathfrak{x}\right),\hfill \\ \hfill \mathbf{W}(\mathfrak{x},0)& ={\mathfrak{c}}_0+{\mathfrak{c}}_{2}tan{h}^2\left(\alpha \mathfrak{x}\right),\hfill \end{array}$$

(46)

where α, $\mathfrak{p}$, ${\mathfrak{c}}_0$ and ${\mathfrak{c}}_2$ are arbitrary constants. For β = 1, the exact solution of (45) is

$$\begin{array}{cc}\hfill \mathbf{U}(\mathfrak{x},\mathfrak{t})& =\frac{\mathfrak{p}-8{\alpha }^2}{3}+4{\alpha }^{2}tan{h}^2(\alpha \mathfrak{x}+\alpha \mathfrak{p}\mathfrak{t}),\hfill \\ \hfill \mathbf{V}(\mathfrak{x},\mathfrak{t})& =-\frac{4{\alpha }^2(3{\alpha }^2{\mathfrak{c}}_0-2\mathfrak{p}{\mathfrak{c}}_2+4{\alpha }^2{\mathfrak{c}}_2)}{3{\mathfrak{c}}_2^2}+\frac{4{\alpha }^4}{{\mathfrak{c}}_2}tan{h}^2(\alpha \mathfrak{x}+\alpha \mathfrak{p}\mathfrak{t}),\hfill \\ \hfill \mathbf{W}(\mathfrak{x},\mathfrak{t})& ={\mathfrak{c}}_0+{\mathfrak{c}}_{2}tan{h}^2(\alpha \mathfrak{x}+\alpha \mathfrak{p}\mathfrak{t}).\hfill \end{array}$$

(47)

Solution: Using the basic theory of the proposed method mentioned in Section 3, the following are the different order problems for the current system.

First-order problem:

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_1(\mathfrak{x},\mathfrak{t})\right\}+{\mathbf{U}}_0(\mathfrak{x},\mathfrak{t})-\left(\frac{1}{s}\right)\left(\frac{\mathfrak{p}-8{\alpha }^2}{3}+4{\alpha }^{2}tan{h}^2\left(\alpha \mathfrak{x}\right)\right)-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{\frac{1}{2}\frac{{\partial }^3{\mathbf{U}}_0}{\partial {\mathfrak{x}}^3}-\\ \hfill 3{\mathbf{U}}_0\frac{\partial {\mathbf{U}}_0}{\partial \mathfrak{x}}+3\frac{\partial \left({\mathbf{V}}_0{\mathbf{W}}_0\right)}{\partial \mathfrak{x}}\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_1(\mathfrak{x},\mathfrak{t})\right\}+{\mathbf{V}}_0(\mathfrak{x},\mathfrak{t})-\left(\frac{1}{s}\right)\left(-\frac{4{\alpha }^2(3{\alpha }^2{\mathfrak{c}}_0-2\mathfrak{p}{\mathfrak{c}}_2+4{\alpha }^2{\mathfrak{c}}_2)}{3{\mathfrak{c}}_2^2}+\frac{4{\alpha }^4}{{\mathfrak{c}}_2}tan{h}^2\left(\alpha \mathfrak{x}\right)\right)-\\ \hfill \left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{{\partial }^3{\mathbf{V}}_0}{\partial {\mathfrak{x}}^3}+3{\mathbf{U}}_0\frac{\partial {\mathbf{V}}_0}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{W}}_1(\mathfrak{x},\mathfrak{t})\right\}+{\mathbf{V}}_0(\mathfrak{x},\mathfrak{t})-\left(\frac{1}{s}\right)\left({\mathfrak{c}}_0+{\mathfrak{c}}_{2}tan{h}^2\left(\alpha \mathfrak{x}\right)\right)-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\{-\frac{{\partial }^3{\mathbf{W}}_0}{\partial {\mathfrak{x}}^3}+\\ \hfill 3{\mathbf{U}}_0\frac{\partial {\mathbf{W}}_0}{\partial \mathfrak{x}}\}=0,\\ \hfill {\mathbf{U}}_1(\mathfrak{x},0)=0,{\mathbf{V}}_1(\mathfrak{x},0)=0,{\mathbf{W}}_1(\mathfrak{x},0)=0.\end{array}$$

(48)

Second-order problem:

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_2(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{\frac{1}{2}\frac{{\partial }^3{\mathbf{U}}_1}{\partial {\mathfrak{x}}^3}-3{\mathbf{U}}_1\frac{\partial {\mathbf{U}}_1}{\partial \mathfrak{x}}+3\frac{\partial \left({\mathbf{V}}_1{\mathbf{W}}_1\right)}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_2(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{{\partial }^3{\mathbf{V}}_1}{\partial {\mathfrak{x}}^3}+3{\mathbf{U}}_1\frac{\partial {\mathbf{V}}_1}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{W}}_2(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{{\partial }^3{\mathbf{W}}_1}{\partial {\mathfrak{x}}^3}+3{\mathbf{U}}_1\frac{\partial {\mathbf{W}}_1}{\partial \mathfrak{x}}\right\}=0,\\ \hfill {\mathbf{U}}_2(\mathfrak{x},0)=0,{\mathbf{V}}_2(\mathfrak{x},0)=0,{\mathbf{W}}_2(\mathfrak{x},0)=0.\end{array}$$

(49)

Third-order problem:

$$\begin{array}{c}\hfill \mathcal{L}\left\{{\mathbf{U}}_3(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{\frac{1}{2}\frac{{\partial }^3{\mathbf{U}}_2}{\partial {\mathfrak{x}}^3}-3{\mathbf{U}}_2\frac{\partial {\mathbf{U}}_2}{\partial \mathfrak{x}}+3\frac{\partial \left({\mathbf{V}}_2{\mathbf{W}}_2\right)}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{V}}_3(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{{\partial }^3{\mathbf{V}}_2}{\partial {\mathfrak{x}}^3}+3{\mathbf{U}}_2\frac{\partial {\mathbf{V}}_2}{\partial \mathfrak{x}}\right\}=0,\\ \hfill \mathcal{L}\left\{{\mathbf{W}}_3(\mathfrak{x},\mathfrak{t})\right\}-\left(\frac{1}{s^{\beta }}\right)\mathcal{L}\left\{-\frac{{\partial }^3{\mathbf{W}}_2}{\partial {\mathfrak{x}}^3}+3{\mathbf{U}}_2\frac{\partial {\mathbf{W}}_2}{\partial \mathfrak{x}}\right\}=0,\\ \hfill {\mathbf{U}}_3(\mathfrak{x},0)=0,{\mathbf{V}}_3(\mathfrak{x},0)=0,{\mathbf{W}}_3(\mathfrak{x},0)=0.\end{array}$$

(50)

Application of inverse Laplace transform gives the solution components of current problems. The fifth-order approximate solution in this case will be represented as

$$\begin{array}{c}\hfill \tilde{\mathbf{U}}=\sum _{m=0}^5{\mathbf{U}}_m(\mathfrak{x},\mathfrak{t}),\\ \hfill \tilde{\mathbf{V}}=\sum _{m=0}^5{\mathbf{V}}_m(\mathfrak{x},\mathfrak{t}),\\ \hfill \tilde{\mathbf{W}}=\sum _{m=0}^5{\mathbf{W}}_m(\mathfrak{x},\mathfrak{t}).\end{array}$$

(51)

while the residual errors of system (45) are

$$\begin{array}{cc}\hfill \mathfrak{Res}\mathbf{U}& =\frac{1}{2}\frac{{\partial }^3\tilde{\mathbf{U}}}{\partial {\mathfrak{x}}^3}-3\tilde{\mathbf{U}}\frac{\partial \tilde{\mathbf{U}}}{\partial \mathfrak{x}}+3\frac{\partial \left(\tilde{\mathbf{V}}\tilde{\mathbf{W}}\right)}{\partial \mathfrak{x}},\hfill \\ \hfill \mathfrak{Res}\mathbf{V}& =-\frac{{\partial }^3\tilde{\mathbf{V}}}{\partial {\mathfrak{x}}^3}+3\tilde{\mathbf{U}}\frac{\partial \tilde{\mathbf{V}}}{\partial \mathfrak{x}},\hfill \\ \hfill \mathfrak{Res}\mathbf{W}& =-\frac{{\partial }^3\tilde{\mathbf{W}}}{\partial {\mathfrak{x}}^3}+3\tilde{\mathbf{U}}\frac{\partial \tilde{\mathbf{W}}}{\partial \mathfrak{x}}.\hfill \end{array}$$

(52)

5. Discussion

This paper is based on the solution and analysis of nonlinear coupled time-fractional KdV systems through a mixed method (Laplace trasform is hybridized with a homotopy perturbation algorithm). The proposed method is applied to different fractional KdV systems (namely a coupled nonlinear fractional KdV system, dispersive long wave time-fractional coupled KdV system, and time-fractional generalized Hiota–Satsuma coupled KdV system), and the obtained results are compared with the already existing literature. Firstly, the proposed method is applied to Example 1, and the obtained results are compared with LADM when the fractional parameter $\beta =1$ in

Table 1. Comparison of LADM and LHPM errors for Example 1 when $\beta$ = 1, $\alpha$ = 1 and $\eta$ = $\mathfrak{a}$ = 0.5. Here, ${\mathbf{E}}_u$ and ${\mathbf{E}}_v$ represents errors of $\mathbf{U}$ and $\mathbf{V}$, respectively.

		Exact Solution		LHPM Solution		LHPM Error		LADM Error [39]
$\mathfrak{x}$	$\mathfrak{t}$	$\mathbf{U}$	$\mathbf{V}$	$\mathbf{U}$	$\mathbf{V}$	${\mathbf{E}}_u$	${\mathbf{E}}_v$	${\mathbf{E}}_u$	${\mathbf{E}}_v$
	0.1	0.01755	0.00877	0.01755	0.00877	5.69 $\times {10}^{-13}$	2.84 $\times {10}^{-13}$	7.44 $\times {10}^{-6}$	4.92 $\times {10}^{-7}$
−10	0.3	0.01734	0.00867	0.01734	0.00867	4.61 $\times {10}^{-11}$	2.30 $\times {10}^{-11}$	2.23 $\times {10}^{-5}$	1.47 $\times {10}^{-6}$
	0.5	0.01713	0.00856	0.01713	0.00856	3.55 $\times {10}^{-10}$	1.77 $\times {10}^{-10}$	3.72 $\times {10}^{-5}$	2.46 $\times {10}^{-6}$
	0.1	0.19717	0.09858	0.19717	0.09858	3.22 $\times {10}^{-12}$	1.61 $\times {10}^{-12}$	3.97 $\times {10}^{-5}$	2.62 $\times {10}^{-6}$
0	0.3	0.19830	0.09915	0.19830	0.09915	2.57 $\times {10}^{-10}$	1.28 $\times {10}^{-10}$	1.19 $\times {10}^{-4}$	7.88 $\times {10}^{-6}$
	0.5	0.19943	0.09971	0.19943	0.09971	1.95 $\times {10}^{-9}$	9.78 $\times {10}^{-10}$	1.98 $\times {10}^{-4}$	1.31 $\times {10}^{-5}$
	0.1	0.00248	0.00124	0.00248	0.00124	1.45 $\times {10}^{-13}$	9.27 $\times {10}^{-14}$	1.07 $\times {10}^{-6}$	7.10 $\times {10}^{-8}$
10	0.3	0.00251	0.00125	0.00251	0.00125	1.18 $\times {10}^{-11}$	5.90 $\times {10}^{-12}$	3.22 $\times {10}^{-6}$	2.13 $\times {10}^{-7}$
	0.5	0.00254	0.00127	0.00254	0.00127	9.13 $\times {10}^{-11}$	4.56 $\times {10}^{-11}$	5.36 $\times {10}^{-6}$	3.55 $\times {10}^{-7}$

. These results clearly indicate that LHPM is more accurate than LADM.

Table 2. LHPM solutions for different $\beta$ in Example 1 when $\mathfrak{x}$ = 7, $\alpha$ = 1.5, $\eta$ = 0.5 and $\mathfrak{a}$ = 0.7 (where ${\mathbf{R}}_u$, ${\mathbf{R}}_v$ and $\mathbf{R}$ represents residual errors of $\mathbf{U}$, $\mathbf{V}$ and system error, respectively).

$\mathit{\beta }$	$\mathfrak{t}$	U	V	${\mathbf{R}}_{\mathit{u}}$	${\mathbf{R}}_{\mathit{v}}$	R
	0.1	0.006960	0.004118	7.29 $\times {10}^{-10}$	4.31 $\times {10}^{-10}$	5.80 $\times {10}^{-10}$
0.33	0.4	0.007153	0.004232	7.18 $\times {10}^{-9}$	4.24 $\times {10}^{-9}$	5.71 $\times {10}^{-9}$
	0.7	0.007264	0.004297	1.80 $\times {10}^{-8}$	1.06 $\times {10}^{-8}$	1.43  $\times {10}^{-8}$
	1.0	0.007348	0.004347	3.25 $\times {10}^{-8}$	1.92 $\times {10}^{-8}$	2.58 $\times {10}^{-8}$
	0.1	0.006789	0.004016	1.62 $\times {10}^{-12}$	9.62 $\times {10}^{-13}$	1.29 $\times {10}^{-12}$
0.66	0.4	0.007008	0.004146	1.57 $\times {10}^{-10}$	9.32 $\times {10}^{-11}$	1.25 $\times {10}^{-10}$
	0.7	0.007177	0.004246	9.98 $\times {10}^{-10}$	5.90 $\times {10}^{-10}$	7.94 $\times {10}^{-10}$
	1.0	0.007327	0.004335	3.23 $\times {10}^{-9}$	1.91 $\times {10}^{-9}$	2.57 $\times {10}^{-9}$
	0.1	0.006727	0.003980	9.96 $\times {10}^{-15}$	5.89 $\times {10}^{-15}$	7.93 $\times {10}^{-15}$
0.88	0.4	0.006922	0.004095	4.44 $\times {10}^{-12}$	2.63 $\times {10}^{-12}$	3.53 $\times {10}^{-12}$
	0.7	0.007103	0.004202	5.22 $\times {10}^{-11}$	3.09 $\times {10}^{-11}$	4.17 $\times {10}^{-11}$
	1.0	0.007280	0.004306	2.51 $\times {10}^{-10}$	1.48 $\times {10}^{-10}$	1.99 $\times {10}^{-10}$

presents LHPM solutions and errors at different values of fractional parameter $\beta$. These results confirm the convergence of LHPM over the entire fractional domain $(0,1)$ and indicate that the proposed method is efficient and reliable. A graphical analysis of Example 1 is also performed in Figure 1, Figure 2 and Figure 3. Figure 1 and Figure 2 depict the LHPM solution and corresponding error as a 3D image, while Figure 3 shows the effect of fractional parameter $\beta$ on the solution. It is observed that $\beta$ has an indirect (inverse) relationship with $\mathbf{U}$ and $\mathbf{V}$. It is also observed that both of the solution components $\mathbf{U}$ and $\mathbf{V}$ increase with time.

Next, we apply LHPM to a nonlinear time-fractional coupled dispersive long-wave KdV system in Example 2, and the numerical results are shown in

Table 3. Comparison of LHPM solution and error in Example 2 when $\beta$ = 1 and $\eta$ = $\mathfrak{a}$ = 0.5 (where ${\mathbf{E}}_u$ and ${\mathbf{E}}_v$ represents errors of $\mathbf{U}$ and $\mathbf{V}$, respectively).

		Exact Solution		LHPM Solution		LHPM Error		LADM Error [39]
$\mathfrak{x}$	$\mathfrak{t}$	$\mathbf{U}$	$\mathbf{V}$	$\mathbf{U}$	$\mathbf{V}$	${\mathbf{E}}_u$	${\mathbf{E}}_v$	${\mathbf{E}}_u$	${\mathbf{E}}_v$
	0.1	0.01071	−0.99469	0.01071	−0.99469	1.49 $\times {10}^{-10}$	6.06 $\times {10}^{-11}$	1.21 $\times {10}^{-7}$	1.31 $\times {10}^{-4}$
−10	0.3	0.01020	−0.99495	0.01020	−0.99495	1.20 $\times {10}^{-8}$	4.88 $\times {10}^{-9}$	6.35 $\times {10}^{-7}$	3.84 $\times {10}^{-4}$
	0.5	0.00970	−0.99519	0.00970	−0.99519	9.21 $\times {10}^{-8}$	3.75 $\times {10}^{-8}$	5.30 $\times {10}^{-7}$	6.27 $\times {10}^{-4}$
	0.1	0.61656	−0.88179	0.61656	−0.88179	1.69 $\times {10}^{-9}$	1.11 $\times {10}^{-9}$	1.23 $\times {10}^{-5}$	7.04 $\times {10}^{-4}$
0	0.3	0.60467	−0.88047	0.60467	−0.88047	1.35 $\times {10}^{-7}$	9.32 $\times {10}^{-8}$	1.14 $\times {10}^{-4}$	2.01 $\times {10}^{-3}$
	0.5	0.59266	−0.87929	0.59266	−0.87929	1.02 $\times {10}^{-6}$	7.38 $\times {10}^{-7}$	3.26 $\times {10}^{-4}$	3.18 $\times {10}^{-3}$
	0.1	0.99582	−0.99792	0.99582	−0.99792	6.25 $\times {10}^{-11}$	2.91 $\times {10}^{-11}$	2.50 $\times {10}^{-6}$	5.08 $\times {10}^{-5}$
10	0.3	0.99561	−0.99781	0.99561	−0.99781	5.11 $\times {10}^{-9}$	2.38 $\times {10}^{-9}$	2.27 $\times {10}^{-5}$	1.56 $\times {10}^{-4}$
	0.5	0.99539	−0.99770	0.99539	−0.99770	3.98 $\times {10}^{-8}$	1.85 $\times {10}^{-8}$	6.36 $\times {10}^{-5}$	2.68 $\times {10}^{-4}$

and

Table 4. LHPM solutions and errors at different values of $\beta$ in Example 2 when $\mathfrak{x}$ = 8 and $\eta$ = $\mathfrak{a}$ = 0.5 (where ${\mathbf{R}}_u$ and ${\mathbf{R}}_v$ represent residual errors in first and second equations, and $\mathbf{R}$ is the system error).

$\mathit{\beta }$	$\mathfrak{t}$	U	V	${\mathbf{R}}_{\mathit{u}}$	${\mathbf{R}}_{\mathit{v}}$	R
	0.1	0.987411	−0.993786	3.33 $\times {10}^{-8}$	1.31 $\times {10}^{-8}$	2.32 $\times {10}^{-8}$
0.33	0.4	0.986286	−0.993240	3.25 $\times {10}^{-7}$	1.26 $\times {10}^{-7}$	2.25 $\times {10}^{-7}$
	0.7	0.985589	−0.992903	8.12 $\times {10}^{-7}$	3.13 $\times {10}^{-7}$	5.62 $\times {10}^{-7}$
	1.0	0.985033	−0.992635	1.45 $\times {10}^{-6}$	5.57 $\times {10}^{-7}$	1.00 $\times {10}^{-6}$
	0.1	0.988326	−0.994231	1.06 $\times {10}^{-10}$	3.18 $\times {10}^{-11}$	6.92 $\times {10}^{-11}$
0.66	0.4	0.987185	−0.993675	1.01 $\times {10}^{-8}$	2.93 $\times {10}^{-9}$	6.54 $\times {10}^{-9}$
	0.7	0.986236	−0.993214	6.34 $\times {10}^{-8}$	1.78 $\times {10}^{-8}$	4.06 $\times {10}^{-8}$
	1.0	0.985345	−0.992782	2.02 $\times {10}^{-7}$	5.54 $\times {10}^{-8}$	1.29 $\times {10}^{-7}$
	0.1	0.988631	−0.994380	1.25 $\times {10}^{-12}$	1.12 $\times {10}^{-13}$	6.84 $\times {10}^{-13}$
0.88	0.4	0.987658	−0.993905	5.46 $\times {10}^{-10}$	3.90 $\times {10}^{-11}$	2.93 $\times {10}^{-10}$
	0.7	0.986699	−0.993438	6.26 $\times {10}^{-9}$	3.36 $\times {10}^{-10}$	3.30 $\times {10}^{-9}$
	1.0	0.985715	−0.992960	2.93 $\times {10}^{-8}$	1.04 $\times {10}^{-9}$	1.52 $\times {10}^{-8}$

. Table 3 presents the comparison of LHPM and LADM at $\beta$ = 1. Analysis of this table showed the dominance of LHPM over the LADM. For the convergence confirmation of the proposed method over the entire fractional domain $(0,1)$, Table 4 displays the solutions along with residual errors at different values of $\beta$. This table indicates that the proposed method is consistent throughout the fractional domain. A graphical analysis of Example 2 is depicted in Figure 4 and Figure 5. These figures show the LHPM solution and error, respectively. The effect of the fractional parameter on the solution components can be seen in Figure 6. It is seen that $\beta$ has a direct relation with solution component $\mathbf{U}$, while it has an indirect relationship with $\mathbf{V}$. It is also noted that $\mathbf{U}$ decreases while $\mathbf{V}$ increases with time.

In the next phase of analysis, a time-fractional generalized Hirota–Satsuma coupled KdV system is computed through LHPM, and a comparison of LHPM and RDTM at $\beta$ = 1 has been made for validity purposes in

Table 5. Comparison of LHPM and RDTM errors in Example 3 when $\beta$ = 1, $\alpha$ = 0.1, ${\mathfrak{c}}_0$ = 1.5 and ${\mathfrak{c}}_2$ = $\mathfrak{p}$ = 0.5 (where ${\mathbf{E}}_u$, ${\mathbf{E}}_v$ and ${\mathbf{E}}_w$ represent errors of $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{W}$, respectively).

		RDTM Error [40]			LHPM Error
$\mathfrak{x}$	$\mathfrak{t}$	${\mathbf{E}}_u$	${\mathbf{E}}_v$	${\mathbf{E}}_w$	${\mathbf{E}}_u$	${\mathbf{E}}_v$	${\mathbf{E}}_w$
	0.1	1.66 $\times {10}^{-11}$	3.32 $\times {10}^{-13}$	2.07 $\times {10}^{-10}$	2.34 $\times {10}^{-16}$	4.69 $\times {10}^{-18}$	2.93 $\times {10}^{-15}$
0.2	0.4	4.24 $\times {10}^{-9}$	8.49 $\times {10}^{-11}$	5.30 $\times {10}^{-8}$	9.59 $\times {10}^{-13}$	1.91 $\times {10}^{-14}$	1.19 $\times {10}^{-11}$
	0.7	3.97 $\times {10}^{-8}$	7.95 $\times {10}^{-10}$	4.96 $\times {10}^{-7}$	2.75 $\times {10}^{-11}$	5.50 $\times {10}^{-13}$	3.43 $\times {10}^{-10}$
	1.0	1.65 $\times {10}^{-7}$	3.30 $\times {10}^{-9}$	2.06 $\times {10}^{-6}$	2.33 $\times {10}^{-10}$	4.66 $\times {10}^{-12}$	2.91 $\times {10}^{-9}$
	0.1	1.63 $\times {10}^{-11}$	3.26 $\times {10}^{-13}$	2.03 $\times {10}^{-10}$	2.27 $\times {10}^{-16}$	4.54 $\times {10}^{-18}$	2.84 $\times {10}^{-15}$
0.5	0.4	4.16 $\times {10}^{-9}$	8.32 $\times {10}^{-11}$	5.20 $\times {10}^{-8}$	9.28 $\times {10}^{-13}$	1.85 $\times {10}^{-14}$	1.16 $\times {10}^{-11}$
	0.7	3.89 $\times {10}^{-8}$	7.78 $\times {10}^{-10}$	4.86 $\times {10}^{-7}$	2.65 $\times {10}^{-11}$	5.31 $\times {10}^{-13}$	3.31 $\times {10}^{-10}$
	1.0	1.61 $\times {10}^{-7}$	3.23 $\times {10}^{-9}$	2.01 $\times {10}^{-6}$	2.24 $\times {10}^{-10}$	4.49 $\times {10}^{-12}$	2.81 $\times {10}^{-9}$
	0.1	1.57 $\times {10}^{-11}$	3.15 $\times {10}^{-13}$	1.96 $\times {10}^{-10}$	2.14 $\times {10}^{-16}$	4.28 $\times {10}^{-18}$	2.67 $\times {10}^{-15}$
0.8	0.4	4.01 $\times {10}^{-9}$	8.02 $\times {10}^{-11}$	5.01 $\times {10}^{-8}$	8.72 $\times {10}^{-13}$	1.74 $\times {10}^{-14}$	1.09 $\times {10}^{-11}$
	0.7	3.74 $\times {10}^{-8}$	7.49 $\times {10}^{-10}$	4.68 $\times {10}^{-7}$	2.49 $\times {10}^{-11}$	4.98 $\times {10}^{-13}$	3.11 $\times {10}^{-10}$
	1.0	1.55 $\times {10}^{-7}$	3.10 $\times {10}^{-9}$	1.94 $\times {10}^{-6}$	2.10 $\times {10}^{-10}$	4.21 $\times {10}^{-12}$	2.63 $\times {10}^{-9}$
	0.1	1.52 $\times {10}^{-11}$	3.05 $\times {10}^{-13}$	1.90 $\times {10}^{-10}$	2.02 $\times {10}^{-16}$	4.05 $\times {10}^{-18}$	2.53 $\times {10}^{-15}$
1.0	0.4	3.88 $\times {10}^{-8}$	7.76 $\times {10}^{-11}$	4.85 $\times {10}^{-8}$	8.23 $\times {10}^{-13}$	1.64 $\times {10}^{-14}$	1.02 $\times {10}^{-11}$
	0.7	3.62 $\times {10}^{-8}$	7.24 $\times {10}^{-10}$	4.52 $\times {10}^{-7}$	2.34 $\times {10}^{-11}$	4.69 $\times {10}^{-13}$	2.93 $\times {10}^{-10}$
	1.0	1.49 $\times {10}^{-7}$	2.99 $\times {10}^{-9}$	1.87 $\times {10}^{-6}$	1.98 $\times {10}^{-10}$	3.96 $\times {10}^{-12}$	2.47 $\times {10}^{-9}$

. The obtained results clearly indicate that LHPM is better than RDTM. For convergence of the proposed scheme, the system is solved for different values of the fractional parameter, and solutions along with residual errors are mentioned in

Table 6. LHPM solutions and errors at different values of fractional parameter $\beta$ in Example 3 when $\alpha$ = 0.2, $\mathfrak{x}$ = 5, ${\mathfrak{c}}_0$ = ${\mathfrak{c}}_2$ = 1 and $\mathfrak{p}$ = 0.55 (where ${\mathbf{R}}_u$, ${\mathbf{R}}_v$, ${\mathbf{R}}_w$ and $\mathbf{R}$ represent residual errors of $\mathbf{U}$, $\mathbf{V}$, $\mathbf{W}$ and system error, respectively).

$\mathit{\beta }$	$\mathfrak{t}$	U	V	W	${\mathbf{R}}_{\mathit{u}}$	${\mathbf{R}}_{\mathit{v}}$	${\mathbf{R}}_{\mathit{w}}$	R
	0.1	0.17473	0.04765	1.61294	6.02 $\times {10}^{-7}$	7.04 $\times {10}^{-9}$	1.10 $\times {10}^{-6}$	5.70 $\times {10}^{-7}$
0.36	0.4	0.17784	0.04778	1.63253	7.31 $\times {10}^{-6}$	8.00 $\times {10}^{-8}$	1.25 $\times {10}^{-5}$	6.63 $\times {10}^{-6}$
	0.7	0.17950	0.04784	1.64306	2.00 $\times {10}^{-5}$	2.10 $\times {10}^{-7}$	3.29 $\times {10}^{-5}$	1.77 $\times {10}^{-5}$
	1.0	0.18069	0.04789	1.65068	3.80 $\times {10}^{-5}$	3.89 $\times {10}^{-7}$	6.08 $\times {10}^{-5}$	3.31 $\times {10}^{-5}$
	0.1	0.17215	0.04755	1.59679	2.32 $\times {10}^{-9}$	2.41 $\times {10}^{-11}$	3.78 $\times {10}^{-9}$	2.04 $\times {10}^{-9}$
0.66	0.4	0.17597	0.04770	1.62071	2.25 $\times {10}^{-7}$	2.13 $\times {10}^{-9}$	3.33 $\times {10}^{-7}$	1.86 $\times {10}^{-7}$
	0.7	0.17867	0.04781	1.63762	1.42 $\times {10}^{-6}$	1.25 $\times {10}^{-8}$	1.95 $\times {10}^{-6}$	1.13 $\times {10}^{-6}$
	1.0	0.18088	0.04790	1.65152	4.64 $\times {10}^{-6}$	3.79 $\times {10}^{-8}$	5.92 $\times {10}^{-6}$	3.53 $\times {10}^{-6}$
	0.1	0.17101	0.04750	1.58966	1.33 $\times {10}^{-11}$	1.86 $\times {10}^{-13}$	2.91 $\times {10}^{-11}$	1.42 $\times {10}^{-11}$
0.88	0.4	0.17458	0.04765	1.61198	5.96 $\times {10}^{-9}$	7.68 $\times {10}^{-11}$	1.20 $\times {10}^{-8}$	6.01 $\times {10}^{-9}$
	0.7	0.17767	0.04777	1.63132	7.00 $\times {10}^{-8}$	8.32 $\times {10}^{-10}$	1.30 $\times {10}^{-7}$	6.70 $\times {10}^{-8}$
	1.0	0.18048	0.04788	1.64892	3.36 $\times {10}^{-7}$	3.68 $\times {10}^{-9}$	5.75 $\times {10}^{-7}$	3.05 $\times {10}^{-7}$

. For graphical analysis, 3D plots of the LHPM solution (see Figure 7) and their error (see Figure 8) are presented to check the behavior of solution components. It is observed that the fractional parameter has an indirect relation with all the solution components (see Figure 9). In addition, $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{W}$ increase with the increase in time.

In all three examples, it is observed that error is comparatively larger near the turbulence/disturbance in the waves. In addition, near the origin of turbulence, the system is more unpredictable, but with time, when waves disperse from the origin, it becomes uniform and the solution is more predictable (see Figure 2, Figure 5 and Figure 8). In addition, the fractional parameter has shown an indirect relationship with the solution components in Examples 1 and 3 (see Figure 3 and Figure 9), but it showed different behavior in Example 2 where the fractional parameter has a direct relationship with $\mathbf{U}$ and an indirect relationship with $\mathbf{V}$ (see Figure 6). Analysis of these figures showed that the behavior of the fractional parameters is not similar and can vary with respect to different KdV systems and conditions.

6. Conclusions

This article is focused on the hybridization of Laplace transform with a homotopy perturbation algorithm for analyzing coupled nonlinear time-fractional KdV systems. Time derivative is considered to be fractional in the Caputo sense in this manuscript. The proposed algorithm is applied to different coupled KdV systems such as dispersive long-wave and generalized Hirota–Satsuma systems. For validity purposes, the obtained solutions are compared with existing ones from the literature, and they concluded that the proposed method is better in terms of accuracy as compared to LADM and RDTM in time-fractional KdV systems. The convergence of the proposed method on the entire fractional domain is confirmed by finding solutions and errors at different values of fractional parameters in the interval $(0,1)$. Numerical and graphical analysis clearly affirm the potential of the proposed technique for dealing with highly nonlinear complex fractional systems. Hence, it is recommended to utilize this approach with Caputo–Fabrizio fractional derivatives and Atangana–Baleanu fractional derivatives having multiple kernels (singular and non-singular) for various classes of higher-order two and three-soliton fractional systems to obtain better results in terms of accuracy with fewer computations in different phenomena of physics and engineering.