1 Introduction

Over the last decades, Fractional Partial Differential Equations (FPDEs) have been a useful tool due to their wide uses for describing the natural phenomena of science and engineering models where the nonlinear wave theory frequently explores asymptotic regimes (such as varying over many scales, high frequency or large amplitude) which are not easily accessible via numerical simulations (see, e.g., Dehghan and Abbaszadeh 2016; Mohebbi Ghandehari and Ranjbar 2013; Thabet and Kendre 2018a, b; Thabet et al 2019b). Solitary waves were discovered by the naval architect John Scott Russell in 1834. When a canal barge hit an underwater obstruction and stopped suddenly, Russell expected that the bow wave would dissolve into lots of little ripples through dispersion. Instead, a smooth, bell-shaped crest perhaps half a meter tall, independent of the cross-channel direction, emerged from the froth. The system of FNPDEs have been increasingly used to represent physical and control systems (see, e.g., Magin et al 2011; Mamchuev 2012, 2016; Thabet et al 2019a; Thabet and Kendre 2019 and references cited therein). Bakkyaraj and Sahadevan (2015) applied the homotopy analysis method to obtain the approximate analytical solution of two coupled time-fractional nonlinear Schrödinger equations. Ouyang (2011) obtained some conditions for the existence of the solutions of a class of nonlinear fractional partial differential equations with delay. Ke et al (2015) proposed the fast direct method for solving the linear block lower triangular Toeplitz-like with tridiagonal blocks system which arises from the time-fractional partial differential equation. Singla and Gupta (2017) introduced an extension of the concept of nonlinear self-adjointness and Noether operators for calculating conserved vectors of the time-fractional nonlinear systems of partial differential equations. Recently, Thabet et al (2017) proposed a new analytical technique for solving a system of fractional nonlinear partial differential equations. Since the exact solutions to a large class of fractional partial differential equations are rarely available, so approximate and numerical methods are applicable. Therefore, accurate methods for finding the solutions of FPDEs are yet under investigation.

The objective of this study is to demonstrate that a general system of FNPDEs can be solved easily by a new approximate-analytical approach and it gives good results in analytical and numerical experiments. The rest of the paper is organized as follows. In Sect. 2, we present basic definitions and preliminaries which are needed in the sequel. In Sect. 3, we introduce a new approximate-analytical approach for solving general systems of FNPDEs. Solitary wave solutions and traveling wave solutions for the systems of time-fractional nonlinear wave equations are obtained in Sect. 4. In Sect. 5, we use Mathematica and Maple software to obtain the numerical solutions in forms of tables and graphs.

2 Basic definitions and preliminaries

There are various definitions and theorems of fractional integrals and derivatives. In this section, we give some of these definitions and theorems which are needed in this paper and can be found in Diethelm (2010), Miller and Ross (1993) and Podlubny (1998) and some references cited therein.

Definition 1

A real valued function u(xt), \(x,t\in {\mathbb {R}},t>0\), is said to be in the space \(C_{\mu }\), \(\mu \in {\mathbb {R}}\) if there exists a real number \(p(>\mu )\), such that \(u(x,t)=t^pu_{1}(x,t)\), where \(u_{1}(x,t)\in C({\mathbb {R}}\times [0,\infty ))\), and it is said to be in the space \(C_{\mu }^{m}\), if and only if \(\frac{{\partial }^{m}u(x,t)}{\partial t^m}\in C_{\mu },m\in {\mathbb {N}}\).

Definition 2

Let \(q, t \in {\mathbb {R}}, m-1 \le q <m \in {\mathbb {N}}\), then the Riemann-Liouville time-fractional partial integral of order q for u(xt) is defined as follows:

$$\begin{aligned} {\mathcal {I}}_t^{q}u(x,t)=\frac{1}{\varGamma (q)}\int _0^t(t-\tau )^{q -1}u(x,\tau ) \mathrm{d}\tau ,~~ t>0. \end{aligned}$$
(1)

Definition 3

Let \(q,t\in {\mathbb {R}}\) and \(m-1<q< m\in {\mathbb {N}}\), then the Caputo time fractional partial derivative of order q of the function u(xt) is defined as follows:

$$\begin{aligned} \left\{ \begin{aligned} {\mathcal {D}}_{t}^{q}u(x,t)&= \int _0^t\frac{(t-\tau )^{m-q -1}}{\varGamma (m-q)}\frac{{\partial }^{m}u(x,\tau )}{\partial {\tau }^m}\mathrm{d}\tau , t>0,\\ {\mathcal {D}}_{t}^{q}u(x,t)&={\mathcal {D}}_{t}^{m}u(x,t)=\frac{\partial ^{m}u(x,t)}{\partial t^{m}},~q=m. \end{aligned} \right. \end{aligned}$$
(2)

Theorem 1

Let \(q,t\in {\mathbb {R}}, m-1<q\le m\in {\mathbb {N}}\) and \(t>0\), then

$$\begin{aligned} \left\{ \begin{aligned} {\mathcal {I}}_{t}^{q}{\mathcal {D}}_{t}^{q}u(x,t)&=u(x,t)-\sum _{k=0}^{m-1}\frac{t^{k}}{k!}\frac{{\partial }^ku(x,0^{+})}{\partial t^{k}},\\ {\mathcal {D}}_{t}^{q}{\mathcal {I}}_{t}^{q}u(x,t)&=u(x,t). \end{aligned} \right. \end{aligned}$$
(3)

Theorem 2

Let \(q,t\in {\mathbb {R}}, m-1<q< m\) and \(m\in {\mathbb {N}}\), then

$$\begin{aligned} \left\{ \begin{aligned} {\mathcal {D}}_{t}^{q}t^p&= \frac{\varGamma (p+1)}{\varGamma (p-q +1)}t^{p-q}={\mathcal {D}}_{t}^{q}t^p,~ p>m-1, ~p\in {\mathbb {R}},\\ {\mathcal {D}}_{t}^{q}t^p&=0,~ p\le m-1, ~p\in {\mathbb {R}}. \end{aligned} \right. \end{aligned}$$
(4)

3 New approximate-analytical approach for solving a system of FNPDEs

In this section, we introduce a new approximate-analytical approach based on homotopy method for solving a general system of FNPDEs with initial values of the following form:

$$\begin{aligned} \left\{ \begin{aligned} {\mathcal {D}}_{t}^{q}u_{i}({\bar{x}},t)&+L_{i}{\bar{u}}+N_{i}{\bar{u}}=f_{i}({\bar{x}},t),~m'-1<q< m'\in {\mathbb {N}},\quad \\ {\mathcal {D}}_{t}^{j}u_{i}({\bar{x}},0)&=f_{ij}({\bar{x}}), ~ j=0,1,2,\ldots ,m'-1,~i=1,2,\ldots ,m, \end{aligned} \right. \end{aligned}$$
(5)

for \(t>0\), where \(L_{i}\) and \(N_{i}\) are linear and nonlinear operators respectively of \({\bar{u}}\) and its partial derivatives which might include other fractional derivatives, and \(f_{i}({\bar{x}},t), f_{ij}({\bar{x}})\) are known analytic functions and \({\mathcal {D}}_{t}^{q}\) is the Caputo time-fractional partial derivative of order q, where we define that

$$\begin{aligned} {\bar{u}}&={\bar{u}}({\bar{x}},t)=(u_{1}({\bar{x}},t),u_{2}({\bar{x}},t),\ldots ,u_{n}({\bar{x}},t)),~{\bar{x}}=(x_{1},x_{2},\ldots ,x_{n})\in {{\mathbb {R}}}^n. \end{aligned}$$
(6)

Theorem 3

For \({\bar{u}}=\sum _{k=0}^{\infty }p^k{\bar{u}}_{k} \) and \({\bar{u}}\) is defined by (6), the nonlinear operators \(N_{i}{\bar{u}}\) are satisfied the following property:

$$\begin{aligned} N_{i}{\bar{u}}= N_{i}\sum _{k=0}^{\infty }p^{k}{\bar{u}}_{k}&=\sum _{n=0}^{\infty }\Bigg [\frac{1}{n!}\frac{{\partial }^n}{\partial {p}^n}\bigg [N_{i}\sum _{k=0}^{n}p^{k}{\bar{u}}_{k}\bigg ]_{p=0}\Bigg ]p^n, \end{aligned}$$
(7)

for \(i=1,2,\ldots ,m.\)

Proof

According to the Maclaurin expansion of \(N_{i}\sum _{k=0}^{\infty }{p}^{k}{\bar{u}}_{k}\) with respect to p, we have

$$\begin{aligned} N_{i}{{\bar{u}}}_{p}&=N_{i}\sum _{k=0}^{\infty }{p}^{k}{\bar{u}}_{k}= [N_{i}\sum _{k=0}^{\infty }{p}^{k}{\bar{u}}_{k}]_{{p}=0}\\&\quad +\bigg [\frac{\partial }{\partial {p}}\bigg [N_{i}\sum _{k=0}^{\infty }{p}^{k}{\bar{u}}_{k}\bigg ]_{{p}=0}\bigg ]{p}+\bigg [\frac{1}{2!}\frac{{\partial }^2}{\partial {p}^2}\bigg [N_{i}\sum _{k=0}^{\infty }{p}^{k}{\bar{u}}_{k}\bigg ]_{{p}=0}\bigg ]{p}^2+\cdots \\&=\sum _{n=0}^{\infty }\bigg [\frac{1}{n!}\frac{{\partial }^n}{\partial {{p}}^n}\bigg [N_{i}\sum _{k=0}^{\infty }{p}^{k}{\bar{u}}_{k}\bigg ]_{{p}=0}\bigg ]{p}^n \\&=\sum _{n=0}^{\infty }\Bigg [\frac{1}{n!}\frac{{\partial }^n}{\partial {{p}}^n}\bigg [N_{i}\bigg (\sum _{k=0}^{n}{p}^{k}{\bar{u}}_{k}+\sum _{k=n+1}^{\infty }{p}^{k}{\bar{u}}_{k}\bigg )\bigg ]_{{p}=0}\Bigg ]{p}^n \\&=\sum _{n=0}^{\infty }\Bigg [\frac{1}{n!}\frac{{\partial }^n}{\partial {{p}}^n}\bigg [N_{i}\bigg (\sum _{k=0}^{n}{p}^{k}{\bar{u}}_{k}+{p}^{n+1}{\bar{u}}_{n+1}+{p}^{n+2}{\bar{u}}_{n+2}+\cdots \bigg )\bigg ]_{{p}=0}\Bigg ]{p}^n \\&=\sum _{n=0}^{\infty }\Bigg [\frac{1}{n!}\frac{{\partial }^n}{\partial {{p}}^n}\bigg [N_{i}\sum _{k=0}^{n}{p}^{k}{\bar{u}}_{k}\bigg ]_{{p}=0}\Bigg ]{p}^n\\&\quad +0~\Bigg (\text {since}~\frac{{\partial }^n}{\partial {{p}}^n}\Bigg [N_{i}({p}^{n+1}{\bar{u}}_{n+1}+{p}^{n+2}{\bar{u}}_{n+2}+\cdots )\Bigg ]_{\lambda =0}=0\Bigg ) \\&=\sum _{n=0}^{\infty }\Bigg [\frac{1}{n!}\frac{{\partial }^n}{\partial {{p}}^n}\bigg [N_{i}\sum _{k=0}^{n}{p}^{k}{\bar{u}}_{k}\bigg ]_{{p}=0}\Bigg ]{p}^n,~i=1,2,\ldots ,m. \end{aligned}$$

\(\square \)

Moreover, let the polynomials \(H_{in}=H_{in}(u_{i0},u_{i1},\ldots ,u_{in})\) for \(i=1,2,\ldots m\) be defined as follows:

$$\begin{aligned} H_{in}(u_{i0},u_{i1},\ldots ,u_{in})&=\frac{1}{n!}\frac{{\partial }^{n}}{\partial {p}^{n}}\bigg [N_{i}\sum _{k=0}^{n}{p}^{k}{\bar{u}}_{k}\bigg ]_{p =0}. \end{aligned}$$
(8)

Therefore, by Theorem 3, the nonlinear operators \(N_{i}{\bar{u}}\) for \(i=1, 2,\ldots , m\) have the following expression:

$$\begin{aligned} N_{i}{\bar{u}}= N_{i}\sum _{k=0}^{\infty }p^{k}{\bar{u}}_{k}=\sum _{n=0}^{\infty }H_{in}p^{n}. \end{aligned}$$
(9)

The following theorem introduces a new approximate-analytical approach based on homotopy method in which the analytical solution for the general system of FNPDEs (5) can be easily obtained.

Theorem 4

Let \(m'-1<q<m',~i=1,2,\ldots ,m\), and \( f_{i}({\bar{x}},t), f_{ij}({\bar{x}})\) are known analytic functions. Then the system (5) admits at least a solution given by:

$$\begin{aligned} u_{i}({\bar{x}},t)&=f_{it}^{(-m')}({\bar{x}},t){+}\sum _{{j}=0}^{m'-1}\frac{t^{j}}{{j}!}f_{i{j}}({\bar{x}}){+}\sum _{k=1}^{\infty }\Bigg [\sum _{s=0}^{k} \frac{(-1)^s\varGamma (k{+}1)}{\varGamma (s{+}1)\varGamma (k-s{+}1)} f_{it}^{(sq-(s+1)m')}({\bar{x}},t)\nonumber \\&~~~-\sum _{r=0}^{k-1}\sum _{s=0}^{k-1-r}\frac{(-1)^s\varGamma (k-r)}{\varGamma (s+1)\varGamma (k-r-s)}\bigg ( L_{irt}^{(sq-(s+1)m')} {+}H_{irt}^{(sq-(s+1)m')}\bigg )\Bigg ], \end{aligned}$$
(10)

where \(f_{it}^{(-m')}({\bar{x}},t)\) means the \({m'}\)-times time integral for \(f_{i}({\bar{x}},t)\) and \(L_{irt}^{(sq-(s+1)m')}=\mathcal {I}_{t}^{sq-(s+1)m'}L_{i}{\bar{u}}_{r}\) and \(H_{irt}^{(sq-(s+1)m')}=\mathcal {I}_{t}^{sq-(s+1)m'}H_{i}(u_{i0}, u_{i1}, \ldots , {u}_{ir})\) denote the \((sq-(s+1)m')\)-times time-fractional partial integral for \(L_{i}{\bar{u}}_{r}\) and \(H_{i}(u_{i0}, u_{i1}, \ldots , {u}_{ir})\) respectively.

Proof

We assume that the system (5) has solutions of the following analytic expansions

$$\begin{aligned} u_{i}({\bar{x}},t)&=\sum _{k=0}^{\infty }u_{ik}({\bar{x}},t),~i{=}1,2,\ldots m. \end{aligned}$$
(11)

Next, in view of homotopy technique (He 1999), we consider the following homotopy:

$$\begin{aligned}&{\mathcal {D}}_{t}^{m'}\varPhi _{i} ({\bar{x}},t,p)\nonumber \\ {}&\quad =f_{i}({\bar{x}},t){+}p[{\mathcal {D}}_{t}^{m'}\varPhi _{i} ({\bar{x}},t,p){-}L_{i}{\bar{\varPhi }} ({\bar{x}},t,p){-}N_{i}{\bar{\varPhi }} ({\bar{x}},t,p){-}{\mathcal {D}}_{t}^{q}\varPhi _{i} ({\bar{x}},t,p)], \end{aligned}$$
(12)

where \(p\in [0,1]\) is an embedding parameter, \({\bar{\varPhi }} ({\bar{x}},t,p)=(\varPhi _{1} ({\bar{x}},t,p),\varPhi _{2} ({\bar{x}},t,p), \ldots , \varPhi _{m} ({\bar{x}},t,p))\) and \(\varPhi _{i} ({\bar{x}},t,p),i=1,2,\ldots ,m\), are unknown functions which can be defined as:

$$\begin{aligned} \varPhi _{i} ({\bar{x}},t,p)=u_{i0}({\bar{x}},t)+\sum _{k=1}^{\infty }p^k u_{ik}({\bar{x}},t),~i=1,2,\ldots ,m. \end{aligned}$$
(13)

When \(p=0\), \(\varPhi _{i}({\bar{x}},t,0)=u_{i0}({\bar{x}},t)\) and when \(p =1\), \(\varPhi _{i}({\bar{x}},t,1)=u_{i}({\bar{x}},t)\). Thus, as p increases from 0 to 1, the solutions \(\varPhi _{i} ({\bar{x}},t,p)\) vary from \(u_{i0}({\bar{x}},t)\) to the solutions \(u_{i}({\bar{x}},t)\).

However, by substituting the system (13) in the system (12), we obtain

$$\begin{aligned} {\mathcal {D}}_{t}^{m'}\sum _{k=0}^{\infty }p^{k}u_{ik}&=f_{i}({\bar{x}},t){+}p\bigg [ {\mathcal {D}}_{t}^{m'}\sum _{k=0}^{\infty }p^{k}u_{ik}{-}L_{i}\sum _{k=0}^{\infty }p^{k}{\bar{u}}_{k}{-}N_{i}\sum _{k=0}^{\infty }p^{k}{\bar{u}}_{k} {-}{\mathcal {D}}_{t}^{q}\sum _{k=0}^{\infty }p^{k}u_{ik} \bigg ], \end{aligned}$$
(14)

for \(i=1,2,\ldots ,m\). From Theorem 3 and system (9), the system (14) can be rewritten as follows:

$$\begin{aligned} {\mathcal {D}}_{t}^{m'}\sum _{k=0}^{\infty }p^{k}u_{ik}&=f_{i}({\bar{x}},t)+p\bigg [ {\mathcal {D}}_{t}^{m'}\sum _{k=0}^{\infty }p^{k}u_{ik}-\sum _{k=0}^{\infty }p^{k}L_{i}{\bar{u}}_{k}-\sum _{n=0}^{\infty }p^nH_{in} -{\mathcal {D}}_{t}^{q}\sum _{k=0}^{\infty }p^{k}u_{ik} \bigg ], \end{aligned}$$
(15)

for \(i=1,2,\ldots ,m\). By equating the terms in the system (15) with identical powers of p, we obtain

$$\begin{aligned} {\mathcal {D}}_{t}^{m'}u_{i0}({\bar{x}},t)&=f_{i}({\bar{x}},t),~~{\mathcal {D}}_{t}^{j}u_{i0}({\bar{x}},0)=f_{ij}({\bar{x}}), i = 1,2,\ldots , m, j = 0,1,\ldots ,m'-1, \end{aligned}$$
(16)

and

$$\begin{aligned} \left\{ \begin{aligned} {\mathcal {D}}_{t}^{m'}u_{i1}&={\mathcal {D}}_{t}^{m'}u_{i0}-L_{i}{\bar{u}}_{0}-H_{i0}-{{\mathcal {D}}_{t}^{q}}u_{i0},~{\mathcal {D}}_{t}^{j}u_{i1}({\bar{x}},0)=0, \\ {\mathcal {D}}_{t}^{m'}u_{i2}&={\mathcal {D}}_{t}^{m'}u_{i1}-L_{i}{\bar{u}}_{1}-H_{i1}-{{\mathcal {D}}_{t}^{q}}u_{i1},~{\mathcal {D}}_{t}^{j}u_{i2}({\bar{x}},0)=0, \\ \vdots&\\ {\mathcal {D}}_{t}^{m'}u_{ik}&={\mathcal {D}}_{t}^{m'}u_{i(k-1)}-L_{i}{\bar{u}}_{k-1}-H_{i(k-1)}-{{\mathcal {D}}_{t}^{q}}u_{i(k-1)},~ {\mathcal {D}}_{t}^{j}u_{ik}({\bar{x}},0)=0, \end{aligned} \right. \end{aligned}$$
(17)

for \(k=1,2,\ldots \) and \(i=1,2,\ldots ,m\) where \({\mathcal {D}}_{t}^{j}u_{ik}({\bar{x}},0)=0\) are obtained by substituting (13) in the initial conditions given by (5). For convenience, we define \(L_{i}{\bar{u}}_{k-1}=L_{i(k-1)}\) and by \(m'\)-fold integral for both sides of the systems (17), with respect to t, we obtain

$$\begin{aligned} u_{i0}&=\sum _{j=0}^{m'-1}\frac{t^{j}}{{j}!} f_{ij}({\bar{x}})+f_{it}^{(-m')}({\bar{x}},t), i=1,2,\ldots , m, \end{aligned}$$
(18)

and

$$\begin{aligned} \left\{ \begin{aligned} u_{i1}&={\mathcal {I}}_{t}^{m'}({\mathcal {D}}_{t}^{m'}u_{i0})-L_{i0t}^{(-m')}-H_{i0t}^{(-m')}-{\mathcal {I}}_{t}^{m'}({{\mathcal {D}}_{t}^{q}}u_{i0}),\\ u_{i2}&={\mathcal {I}}_{t}^{m'}{\mathcal {D}}_{t}^{m'}u_{i1}-L_{i1t}^{(-m')}-H_{i1t}^{(-m')}-{\mathcal {I}}_{t}^{m'}({{\mathcal {D}}_{t}^{q}}u_{i1}),\\ \vdots&\\ u_{ik}&={\mathcal {I}}_{t}^{m'}{\mathcal {D}}_{t}^{m'}u_{i(k-1)}-L_{i(k-1)t}^{(-m')}({\bar{x}},t)-H_{i(k-1)t}^{(-m')}-{\mathcal {I}}_{t}^{m'}({{\mathcal {D}}_{t}^{q}}u_{i(k-1)}),~{\mathcal {D}}_{t}^{j}u_{ik}({\bar{x}},0)=0, \end{aligned} \right. \end{aligned}$$
(19)

for \(k=1,2,\ldots ,\) and \(i=1,2,\ldots ,m\). By substituting \(u_{i(k-1)}\) into \(u_{ik}\), \(k=1,2,\ldots \) (for example, substituting \(u_{i0}\) into \(u_{i1}\) to get \(u_{i1}\) and since \({\mathcal {I}}_{t}^{m'}{{\mathcal {D}}_{t}^{q}}u_{i0}\ne {\mathcal {I}}_{t}^{m'-q}u_{i0}\) and by using Theorem 2, we have \( {\mathcal {I}}_{t}^{m'}({{\mathcal {D}}_{t}^{q}}u_{i0})={\mathcal {I}}_{t}^{m'}(\underbrace{\sum _{j=0}^{m'-1}{{\mathcal {D}}_{t}^{q}}\frac{t^{j}}{{j}!} f_{ij}({\bar{x}})}_{0}+{{\mathcal {D}}_{t}^{q}}f_{it}^{(-m')}({\bar{x}},t))={\mathcal {I}}_{t}^{m'}f_{it}^{(q-m')}({\bar{x}},t)=f_{it}^{(q-2m')}({\bar{x}},t)\)). By using Theorem 1 with the given initial conditions, we obtain

$$\begin{aligned} \left\{ \begin{aligned} u_{i1}&=f_{it}^{(-m')}({\bar{x}},t)-L_{i0t}^{(-m')}-H_{i0t}^{(-m')}-f_{it}^{(q-2m')}({\bar{x}},t),\\ u_{i2}&=\sum _{s=0}^{2} \frac{(-1)^s\varGamma (2+1)}{\varGamma (s+1)\varGamma (2-s+1)}f_{it}^{(sq{-}(s+1)m')}({\bar{x}},t){-}\sum _{r=0}^{2-1}\sum _{s=0}^{2-1-r}\frac{(-1)^s\varGamma (2-r)}{\varGamma (s+1)\varGamma (2-r-s)}\\&\quad \times \Bigg [ L_{irt}^{(sq-(s+1)m')}+H_{irt}^{(sq-(s+1)m')}\Bigg ],\\ \vdots&\\ u_{ik}&=\sum _{s=0}^{k} \frac{(-1)^s\varGamma (k+1)}{\varGamma (s+1)\varGamma (k-s+1)}f_{it}^{(sq-(s+1)m')}({\bar{x}},t)\\&\quad -\sum _{r=0}^{k-1}\sum _{s=0}^{k-1-r}\frac{(-1)^s\varGamma (k-r)}{\varGamma (s+1)\varGamma (k-r-s)}\Bigg [ L_{irt}^{(sq-(s+1)m')} +H_{irt}^{(sq-(s+1)m')}\Bigg ], \end{aligned} \right. \end{aligned}$$
(20)

for \(k=1,2,\ldots \) and \(i=1,2,\ldots ,m\). From system (13), we obtain

$$\begin{aligned} u_{i}({\bar{x}},t)=\varPhi _{i} ({\bar{x}},t,1)=u_{i0}({\bar{x}},t)+\sum _{k=1}^{\infty }u_{ik}({\bar{x}},t), \end{aligned}$$
(21)

for \(i=1,2,\ldots m.\) Substituting (18) and the components of (20) in the system (21) completes the proof. \(\square \)

Theorem 5

Let B be a Banach space. Then the series solution given by (11) converges to \(S_{i}\in B\) for \(i=1,2,\ldots ,m,\) if there exists \({\gamma }_{i},~0\le {\gamma }_{i}<1\) such that, \(\Vert u_{in}\Vert \le {\gamma }_{i}\Vert u_{i(n-1)}\Vert \) for \(\forall n\in {\mathbb {N}}\).

Proof

Define the sequences \({ S_{in} }\), \(i=1,2,\ldots , n\) of partial sums of the series given by the system (11) as

$$\begin{aligned}&\left\{ \begin{aligned} S_{i0}&=u_{i0}({\bar{x}},t),\\ S_{i1}&=u_{i0}({\bar{x}},t)+u_{i1}({\bar{x}},t),\\ S_{i2}&=u_{i0}({\bar{x}},t)+u_{i1}({\bar{x}},t)+u_{i2}({\bar{x}},t),\\ \vdots ~~&\\ S_{in}&=u_{i0}({\bar{x}},t)+u_{i1}({\bar{x}},t)+u_{i2}({\bar{x}},t)+\cdots +u_{in}({\bar{x}},t),~i=1,2,\ldots ,m, \end{aligned} \right. \end{aligned}$$
(22)

and we need to show that \(\left\{ S_{in} \right\} \) are Cauchy sequences in a Banach space B. For this purpose, we consider

$$\begin{aligned} \Vert S_{i(n+1)}-S_{in}\Vert&=\Vert u_{i(n+1)}({\bar{x}},t)\Vert {\le } {\gamma }_{i}\Vert u_{in}({\bar{x}},t)\Vert {\le } {\gamma }_{i}^{2}\Vert u_{i(n-1)}({\bar{x}},t)\Vert {\le }{\cdots }{\le } {\gamma }_{i}^{n+1}\Vert u_{i0}({\bar{x}},t)\Vert , \end{aligned}$$
(23)

where \(~i=1,2,\ldots ,m\). For every \(n,n'\in {\mathbb {N}},~n\ge n'\), by using the system (23) and triangle inequality successively, we have,

$$\begin{aligned} \Vert S_{in}-S_{in'}\Vert&=\Vert S_{i(n'+1)}-S_{in'}+S_{i(n'+2)}-S_{i(n'+1)}+\cdots + S_{in}-S_{i(n-1)}\Vert \nonumber \\&\le \Vert S_{i(n'+1)}-S_{in'}\Vert +\Vert S_{i(n'+2)}-S_{i(n'+1)}\Vert +\cdots +\Vert S_{in}-S_{i(n-1)}\Vert \nonumber \\&\le {\gamma }_{i}^{n'+1}\Vert u_{i0}({\bar{x}},t)\Vert +{\gamma }_{i}^{n'+2}\Vert u_{i0}({\bar{x}},t)\Vert +\cdots +{\gamma }_{i}^{n}\Vert u_{i0}({\bar{x}},t)\Vert \nonumber \\&={\gamma }_{i}^{n'+1}(1+{\gamma }_{i}+\cdots +{\gamma }_{i}^{n-n'-1})\Vert u_{i0}({\bar{x}},t)\Vert \nonumber \\&\le {\gamma }_{i}^{n'+1}\Bigg (\frac{1-{\gamma }^{n-n'}}{1-\gamma _{i}}\Bigg )\Vert u_{i0}({\bar{x}},t)\Vert . \end{aligned}$$
(24)

Since \(0<\gamma _{i}<1\), so \(1-{\gamma }_{i}^{n-n'}\le 1\). Then

$$\begin{aligned} \Vert S_{in}-S_{in'}\Vert \le \frac{{\gamma }_{i}^{n'+1}}{1-\gamma }_{i}\Vert u_{i0}({\bar{x}},t)\Vert . \end{aligned}$$
(25)

Since \(u_{i0}({\bar{x}},t)\) is bounded, then

$$\begin{aligned} \lim _{n,n'\rightarrow \infty }\Vert S_{in}-S_{in'}\Vert =0,~i=1,2,\ldots ,m. \end{aligned}$$
(26)

Therefore, the sequences \(\left\{ S_{in} \right\} \) are Cauchy sequences in the Banach space B, so the series of solutions defined by the system (22) converges. This completes the proof. \(\square \)

Theorem 6

The maximum absolute truncation error of the series solution (10) of the nonlinear fractional partial differential system (5) is estimated to be

$$\begin{aligned} \sup _{({\bar{x}},t)\in \varOmega }\bigg |u_{i}({\bar{x}},t)-\sum _{k=0}^{m}u_{ik}({\bar{x}},t)\bigg |\le \frac{{\gamma }_{i}^{m+1}}{1-{\gamma }_{i}}\sup _{({\bar{x}},t)\in \varOmega }|u_{i0}({\bar{x}},t)|, \end{aligned}$$
(27)

for \(i=1,2,\ldots ,m\) where the region \(\varOmega \subset {{\mathbb {R}}}^{n+1}\).

Proof

From Theorem 5, we have

$$\begin{aligned} \Vert S_{in}-S_{in'}\Vert \le \frac{{\gamma }_{i}^{n'+1}}{1-{\gamma }_{i}}\sup _{({\bar{x}},t)\in \varOmega }|u_{i0}({\bar{x}},t)|,~i=1,2,\ldots ,m. \end{aligned}$$
(28)

But we assume that \(S_{in}=\sum _{k=0}^{n}u_{ik}({\bar{x}},t)\) for \(i=1,2,\ldots ,m\), and since \(n\rightarrow \infty \), we obtain \(S_{in}\rightarrow u_{i}({\bar{x}},t)\), so the system (28) can be rewritten as

$$\begin{aligned} \Vert u_{i}({\bar{x}},t)-S_{in'}\Vert&=\Vert u_{i}({\bar{x}},t)-\sum _{k=0}^{n'}u_{ik}({\bar{x}},t)\Vert \le \frac{{\gamma }_{i}^{n'+1}}{1-{\gamma }_{i}}\sup _{({\bar{x}},t)\in \varOmega }|u_{i0}({\bar{x}},t)|,~i=1,2,\ldots ,m. \end{aligned}$$
(29)

So, the maximum absolute truncation error in the region \(\varOmega \) is

$$\begin{aligned} \sup _{({\bar{x}},t)\in \varOmega }\bigg |u_{i}({\bar{x}},t)-\sum _{k=0}^{n'}u_{ik}({\bar{x}},t)\bigg |\le \frac{{\gamma }_{i}^{n'+1}}{1-{\gamma }_{i}}\sup _{({\bar{x}},t)\in \varOmega }|u_{i0}({\bar{x}},t)|,~i=1,2,\ldots ,m. \end{aligned}$$
(30)

and this completes the proof. \(\square \)

4 Solitary and traveling wave solutions for systems of time-fractional nonlinear wave equations

This section includes Solitary and traveling wave solutions for some examples of systems of time-fractional nonlinear wave equations. These examples are chosen because their closed form solution are available, or they have been solved previously by some other well-known methods. We use the following mathematical notations \(u_{x}=\frac{\partial u(x,t)}{\partial x}, u_{xx}=\frac{\partial ^2 u(x,t)}{\partial x2}, u_{xxx}=\frac{\partial ^3 u(x,t)}{\partial x3}\).

Example 1

(Paquin and Winternitz 1990; Thabet et al 2017; Wang 1995) Consider the following system of time-fractional dispersive long wave equations with initial values:

$$\begin{aligned} \left\{ \begin{aligned} {\mathcal {D}}_{t}^{q}u+&v_{x}+uu_{x}=0,~ u(x,0)=\alpha \bigg [1+\tanh \left( \frac{1}{2} [\beta +\alpha x]\right) \bigg ],\\ {\mathcal {D}}_{t}^{q}v+&uv_{x}+vu_{x}+u_{x}+u_{xxx}=0,~v(x,0)=-1+\frac{1}{2} \alpha ^2 \text {sech}^2\left( \frac{1}{2} [\beta +\alpha x]\right) , \end{aligned} \right. \end{aligned}$$
(31)

for \(~0<q<1\). For \(q=1\), the exact solitary and traveling wave solutions for the system (31) are given by

$$\begin{aligned} u_{EX}(x,t)&=\alpha \bigg [1{+}\tanh \left( \frac{1}{2} [\beta {+}\alpha x-{\alpha }^2 t]\right) \bigg ],~ v_{EX}(x,t){=-}1{+}\frac{1}{2} \alpha ^2 \text {sech}^2\left( \frac{1}{2} [\beta {+}\alpha x-{\alpha }^2 t]\right) ,\nonumber \\ \end{aligned}$$
(32)

where \(\alpha ,\beta \) are arbitrary constants.

By comparing the system (31) with the system (5) for \(m'=1\), we obtain

$$\begin{aligned} L_{1}(u,v)= & {} v_{x},~L_{2}(x,t)=u_{x}+u_{xxx},\nonumber \\ N_{1}(u,v)= & {} uu_{x},~N_{2}(u,v)=uv_{x}+vu_{x}, f_{1}(x,t)=0, f_{2}(x,t)=0, \end{aligned}$$
(33)

From (11), we assume that the system (31) has solutions given by

$$\begin{aligned} u(x, t) = \sum _{k=0}^{\infty } u_{k} (x, t), v(x, t) = \sum _{k=0}^{\infty } v_{k} (x, t), \end{aligned}$$
(34)

for \(m'=1\) and \(i=1,2\), we have

$$\begin{aligned} u_{0}(x,t)&=\sum _{j=0}^{1-1}\frac{t^{j}}{{j}!}{\mathcal {D}}_{t}^{j}u_{0}(x,0),~v_{0}(x,t)=\sum _{j=0}^{1-1}\frac{t^{j}}{{j}!}{\mathcal {D}}_{t}^{j}v_{0}(x,0). \end{aligned}$$
(35)

By using (33) with initial conditions given by (31), we obtain

$$\begin{aligned} u_{0}(x,t)&=\alpha \bigg [1+\tanh \left( \frac{1}{2} [\beta +\alpha x]\right) \bigg ],~ v_{0}(x,t)=-1+\frac{1}{2} \alpha ^2 \text {sech}^2\left( \frac{1}{2} [\beta +\alpha x]\right) . \end{aligned}$$
(36)

Next, from (20), we have

$$\begin{aligned} \begin{aligned} u_{k}(x,t)&=-\sum _{r=0}^{k-1}\sum _{s=0}^{k-1-r}\frac{(-1)^s\varGamma (k-r)}{\varGamma (s+1)\varGamma (k-r-s)}\Bigg [ L_{1rt}^{(sq-(s+1))} +H_{1rt}^{(sq-(s+1))}\Bigg ],\\ v_{k}(x,t)&=-\sum _{r=0}^{k-1}\sum _{s=0}^{k-1-r}\frac{(-1)^s\varGamma (k-r)}{\varGamma (s+1)\varGamma (k-r-s)}\Bigg [ L_{2rt}^{(sq-(s+1))} +H_{2rt}^{(sq-(s+1))}\Bigg ]. \end{aligned} \end{aligned}$$
(37)

Consequently, by using Theorem 4, systems (36) and (37) with the help of Mathematica software, the first three components of the wave solutions for the system (31) are derived as follows:

$$\begin{aligned} u_{1}(x,t)= & {} -\frac{\alpha ^3}{2} \text {sech}^2\left( \frac{1}{2} [\beta +\alpha x]\right) t,~ v_{1}(x,t)= 4 \alpha ^4 \sinh ^4\left( \frac{1}{2} [\beta +\alpha x]\right) \text {csch}^3(\beta +\alpha x) t,\nonumber \\ u_{2}(x,t)= & {} \frac{1}{4} \alpha ^3 \text {sech}^2\left( \frac{1}{2} [\beta +\alpha x]\right) \Bigg [-2-\alpha ^2 \tanh \left( \frac{1}{2} [\beta +\alpha x]\right) t+\frac{2 t^{1-q}}{\varGamma (3-q)}\Bigg ]t,\nonumber \\ v_{2}(x,t)= & {} \frac{1}{8} \alpha ^4 \text {sech}^2\left( \frac{1}{2} [\beta +\alpha x]\right) \nonumber \\&\quad \Bigg [-\frac{2 \sinh (\beta +\alpha x) t^{1-q}}{\varGamma (3-q)}+2 \sinh (\beta +\alpha x)+\alpha ^2 [\cosh (\beta +\alpha x)-2]t\Bigg ]t,\nonumber \\ u_{3}(x,t)= & {} \frac{\alpha ^3 (\sinh (\beta +\alpha x)+\cosh (\beta +\alpha x)+1)^6}{96 (\cosh (\beta +\alpha x)+1)^5}\bigg [\cosh (4 [\beta + \alpha x])-\sinh (4 [\beta + \alpha x])\bigg ]\nonumber \\&\quad \times \Bigg [6 t^{1-2 q} \Bigg [\sinh (\beta {+}\alpha x){+}\cosh (\beta {+}\alpha x){+}1\Bigg ]\Bigg [{-}\frac{2 t^q}{\varGamma (4{-}q)} \bigg [[q{-}\alpha ^2 t{-}3] [\sinh (\beta {+}\alpha x)\nonumber \\&\quad +\cosh (\beta +\alpha x)]+q+\alpha ^2 t-3\bigg ]\nonumber \\&\quad -\frac{1}{\varGamma (4-2 q)} \bigg [\sinh (\beta +\alpha x)+\cosh (\beta +\alpha x)+\bigg ]t\Bigg ]-2 \bigg [\sinh (\beta +\alpha x)\nonumber \\&\quad +\cosh (\beta +\alpha x)\bigg ] \bigg [-2 \alpha ^4 t^2+\alpha ^4\nonumber \\&\quad \cosh (\beta +\alpha x)t^2+6 \alpha ^2 \sinh (\beta +\alpha x)t+6 \cosh (\beta +\alpha x)+6\bigg ]\Bigg ]t,\nonumber \\ v_{3}(x,t)= & {} \frac{\alpha ^4 \cosh ^4\left( \frac{1 }{2}[\beta +\alpha x]\right) }{3 (\cosh (\beta +\alpha x)+1)^4}\Bigg [\cosh (\beta +\alpha x)-\sinh (\beta +\alpha x)\Bigg ]\nonumber \\&\quad \bigg [6\bigg [\sinh (2 \beta +2 \alpha x)+\cosh (2 \beta +2 \alpha x)-1\bigg ] \frac{t^{2-2 q} }{\varGamma (4-2 q)}\nonumber \\&\quad +2 \bigg [\sinh (\beta +\alpha x)+\cosh (\beta +\alpha x)\bigg ] \Bigg [\bigg [q \sinh (\beta +\alpha x)\nonumber \\&\quad +2 \alpha ^2[1-\cosh (\beta +\alpha x)]t-3 \sinh (\beta +\alpha x)\bigg ]\frac{12t^{1-q} }{\varGamma (4-q)}\nonumber \\&\quad +\alpha ^4 t^2 \bigg [\cosh (\beta +\alpha x)-5\bigg ] \tanh (\frac{1}{2}[\beta +\alpha x])-12 \alpha ^2 t\nonumber \\&\quad +6 \alpha ^2 \cosh (\beta +\alpha x)t+6 \sinh (\beta +\alpha x)\Bigg ]\bigg ]t,\nonumber \\&\vdots \end{aligned}$$
(38)

and so on. By substituting (36) and (38) into (34), we obtain the third-order approximate wave solutions for (31) as follows:

$$\begin{aligned} u(x,t)&\approx \sum _{k=0}^{3} u_{k}({x},t), ~ v(x,t)\approx \sum _{k=0}^{3} v_{k}({x},t). \end{aligned}$$
(39)

Example 2

(Wang 1995; Wang et al 1996) Consider the following system of time-fractional approximate long water wave equations with initial values:

$$\begin{aligned} \left\{ \begin{aligned} {\mathcal {D}}_{t}^{q}u&-uu_{x}- v_{x}+\frac{1}{2}u_{xx}=0,~ u(x,0)=\frac{1}{2} \bigg [\alpha +2 \beta +\alpha \tanh (\frac{1}{2} [\gamma +\alpha x])\bigg ],\\ {\mathcal {D}}_{t}^{q}v&-(uv)_{x}- \frac{1}{2}v_{xx}=0,~v(x,0)=\frac{1}{4} \alpha ^2 \text {sech}^2\left( \frac{1}{2} \bigg [\gamma +\alpha x\bigg ]\right) , \end{aligned} \right. \end{aligned}$$
(40)

for \(0<q<1\). For \(q=1\), the exact solitary and traveling wave solutions of the system (40) are given by

$$\begin{aligned} u_{EX}(x,t)&=\frac{1}{2} \bigg [\alpha {+}2 \beta {+}\alpha \tanh (\frac{1}{2} [\gamma {+}\alpha x{+}ct])\bigg ],~ v_{EX}(x,t){=}\frac{\alpha ^2}{4}\text {sech}^2(\frac{1}{2} \bigg [\gamma +\alpha x+ct\bigg ]), \end{aligned}$$
(41)

where \(\alpha ,\gamma \) and \(c=(\alpha \beta +\frac{\alpha ^2}{2})\) are arbitrary constants.

Table 1 Numerical values of the wave solutions and absolute error when \(\alpha =\beta =0.50\) and \(q=0.50,0.75,1\) for Example 1
Table 2 Numerical values of the wave solutions and the absolute error when \(\alpha =\beta =0.50, \gamma =1\) and \(q=0.50,0.75,1\) for Example 2

In order to obtain the wave solutions for system (40), we compare the system (40) with the system (5), we get

$$\begin{aligned} L_{1}(u,v)&=\frac{1}{2}u_{xx}-v_{x}, ~L_{2}(u,v)=-\frac{1}{2}v_{xx}, N_{1}(u,v)=-uu_{x},\nonumber \\ N_{2}(u,v)&=-(uv_{x}+vu_{x}), f_{1}(x,t)=0, f_{2}(x,t)=0. \end{aligned}$$
(42)

From (11), we assume that the system (40) has solutions given by

$$\begin{aligned} u(x, t) = \sum _{k=0}^{\infty } u_{k} (x, t), v(x, t) = \sum _{k=0}^{\infty } v_{k} (x, t), \end{aligned}$$
(43)

From (18), for \(m'=1\) and \(i=1,2\), we have

$$\begin{aligned} u_{0}(x,t)&=\sum _{j=0}^{1-1}\frac{t^{j}}{{j}!}{\mathcal {D}}_{t}^{j}u_{0}(x,0),~v_{0}(x,t)=\sum _{j=0}^{1-1}\frac{t^{j}}{{j}!}{\mathcal {D}}_{t}^{j}v_{0}(x,0). \end{aligned}$$
(44)

By using (42) with initial conditions given by (40), we obtain

$$\begin{aligned} u_{0}(x,t)&=\frac{1}{2} \bigg [\alpha +2 \beta +\alpha \tanh (\frac{1}{2} [\gamma +\alpha x])\bigg ],~ v_{0}(x,t)=\frac{1}{4} \alpha ^2 \text {sech}^2(\frac{1}{2} \bigg [\gamma +\alpha x\bigg ]). \end{aligned}$$
(45)

Next, from (20), we have

$$\begin{aligned} \begin{aligned} u_{k}(x,t)&=-\sum _{r=0}^{k-1}\sum _{s=0}^{k-1-r}\frac{(-1)^s\varGamma (k-r)}{\varGamma (s+1)\varGamma (k-r-s)}\Bigg [ L_{1rt}^{(sq-(s+1))} +H_{1rt}^{(sq-(s+1))}\Bigg ],\\ v_{k}(x,t)&=-\sum _{r=0}^{k-1}\sum _{s=0}^{k-1-r}\frac{(-1)^s\varGamma (k-r)}{\varGamma (s+1)\varGamma (k-r-s)}\Bigg [ L_{2rt}^{(sq-(s+1))} +H_{2rt}^{(sq-(s+1))}\Bigg ]. \end{aligned} \end{aligned}$$
(46)

By applying Theorem 4 into using (45) and (46) with the help of Mathematica software, we obtain the first three components of the wave solutions for system (40) as follows:

$$\begin{aligned} \begin{aligned} u_{1}(x,t)&=\frac{\alpha ^2 (\alpha +2 \beta )t}{4 [\cosh (\gamma +\alpha x)+1]},~ v_{1}(x,t)=-\alpha ^3 (\alpha +2 \beta ) \sinh ^4(\frac{1}{2} [\gamma +\alpha x]) \text {csch}^3(\gamma +\alpha x)t,\\ u_{2}(x,t)&=\frac{\alpha ^2 [\alpha +2 \beta ] \cosh ^3(\frac{1}{2} [\gamma +\alpha x])t}{4 (\cosh (\gamma +\alpha x)+1)^3} \bigg [-\frac{1}{\varGamma (3-q)}4 \cosh (\frac{1}{2} [\gamma +\alpha x])t^{1-q}-\alpha (\alpha +2 \beta ) \sinh (\frac{1}{2} [\gamma +\alpha x])t\\&\quad +4 \cosh (\frac{1}{2} [\gamma +\alpha x])\bigg ],\\ v_{2}(x,t)&=\frac{1}{64} \alpha ^3 [\alpha +2 \beta ] \text {sech}^4(\frac{1}{2} [\gamma +\alpha x])t \bigg [\frac{1}{\varGamma (3-q)}4 \sinh (\gamma +\alpha x) t^{1-q}+\alpha [\alpha +2 \beta ] [\cosh (\gamma +\alpha x)\\&\quad -2]t-4 \sinh (\gamma +\alpha x)\bigg ],\\ u_{3}(x,t)&=\frac{1}{384} \alpha ^2 [\alpha +2 \beta ] \text {sech}^4(\frac{1}{2} [\gamma +\alpha x])t \bigg [12 t \Bigg [2 \bigg [\frac{1}{\varGamma (4-q)}\bigg [2 (q-3)\ [\cosh (\gamma +\alpha x)+1]+\alpha t [\alpha +2 \beta ] \\&\quad \times \sinh (\gamma +\alpha x)\bigg ]t^q+\frac{t (\cosh (\gamma +\alpha x)+1)}{\varGamma (4-2 q)}\bigg ]t^{-2 q}-\alpha [\alpha +2 \beta ] \sinh (\gamma +\alpha x)\Bigg ]-2 \bigg [\alpha ^2 t^2 (\alpha +2 \beta )^2\\&\quad -12\bigg ]+\bigg [\alpha ^2 (\alpha +2 \beta )^2 t^2+24\bigg ] \cosh (\gamma +\alpha x)\bigg ],\\ v_{3}(x,t)&=\frac{1}{768} \alpha ^3 [\alpha +2 \beta ] \bigg [-\frac{48 \text {sech}^4(\frac{1}{2} [\gamma +\alpha x])}{\varGamma (4-q)} \Bigg [2 (q-3) \sinh (\gamma +\alpha x)+\alpha [\alpha +2 \beta ] [\cosh (\gamma +\alpha x)\\&\quad -2]t\Bigg ]t^{1-q}-\frac{768 \sinh ^4(\frac{1}{2} [\gamma +\alpha x])t^{2-2 q} }{\varGamma (4-2 q)\text {cosh}^3(\gamma +\alpha x)}+\alpha ^2(\alpha +2 \beta )^2 \Bigg [11 \sinh (\frac{1}{2} [\gamma +\alpha x])-\sinh (\frac{3}{2} [\gamma +\alpha x])\Bigg ]\\&\quad \times \text {sech}^5(\frac{1}{2} [\gamma +\alpha x])t^2+24~ \text {sech}^4(\frac{1}{2} [\gamma +\alpha x]) \bigg [\alpha [\alpha +2 \beta ][\cosh (\gamma +\alpha x)-2]t-2 \sinh (\gamma +\alpha x)\bigg ]\bigg ]t,\\ \vdots&\end{aligned} \end{aligned}$$
(47)

and so on. By substituting (45) and (47) into (43), we obtain the third-order approximate wave solutions for system (40) as follows:

$$\begin{aligned} u(x,t)&\approx \sum _{k=0}^{3} u_{k}({x},t), ~ v(x,t)\approx \sum _{k=0}^{3} v_{k}({x},t). \end{aligned}$$
(48)

5 Discussion and numerical results

In Table 1, we evaluate numerical values of the wave solutions as well as the absolute error for Example 1 among different values of xt when \(\alpha =\beta =0.5\) and \(q=0.50, 0.75, 1\). Table 2 includes numerical values of the wave solutions with absolute error for Example 2 by different values of xt when \(\alpha =\beta =0.5, \gamma =1\) and \(q=0.50, 0.75, 1\).

The approximate wave solutions for Example 1 is graphically represented in Fig. 1 when \(\alpha =\beta =0.50\) and \(q=0.25\), in Fig. 2 when \(\alpha =\beta =0.50\) and \(q=0.50\), in Fig. 3 when \(\alpha =\beta =0.50\) and \(q=0.75\) and in Fig. 4 when \(\alpha =\beta =0.50\) and \(q=1\) respectively among different values of xt. In Fig. 5, we plot the graphs of the exact wave solutions for Example 1. In Fig. 6 we plot the graphs of the wave solutions vs x when t is fixed at \(t=10\) and \(\alpha =\beta =0.50\) among different values of q for Example 1. From Fig. 6a when the fractional order q decreases, the wave starts dissipating accordingly and such types of waves are called travelling waves. From Fig. 6b, we see that the amplitudes of the wave are inversely proportional to the fractional order q and these kinds of waves are called solitary waves. We plot the graphs of the approximate solitary and travelling wave solutions for Example 2 in Fig. 7 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=0.25\), Fig. 8 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=0.50\), Fig. 9 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=0.75\) and in Fig. 10 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=1\) among different values of xt. In Fig. 11, we plot the graphs of the exact wave solutions for Example 2. In Fig. 12 we plot the graphs of the wave solutions vs x when t is fixed at \(t=10\) and \(\alpha =\beta =0.50, \gamma =1\) among different values of q for Example 2. From Fig. 12a, we see that the wave starts dissipating when the fractional order q increases accordingly and such types of waves are called travelling waves. From Fig. 12b, we see that the amplitudes of the wave are inversely proportional to the fractional order q and these kinds of waves are called solitary waves.

Fig. 1
figure 1

The graphs of the approximate wave solutions u(xt), v(xt) for Example 1 when \(\alpha =\beta =0.50\) and \(q=0.25\)

Fig. 2
figure 2

The graphs of the approximate wave solutions u(xt), v(xt) for Example 1 when \(\alpha =\beta =0.50\) and \(q=0.50\)

Fig. 3
figure 3

The graphs of the approximate wave solutions u(xt), v(xt) for Example 1 when \(\alpha =\beta =0.50\) and \(q=0.75\)

Fig. 4
figure 4

The graphs of the approximate wave solution u(xt), v(xt) for Example 1 when \(\alpha =\beta =0.50\) and \(q=1\)

Fig. 5
figure 5

The graphs of the exact solutions u(xt), v(xt) for Example 1 when \(\alpha =\beta =0.50\)

Fig. 6
figure 6

The graphs of the approximate wave solutions u(xt), v(xt) vs. x at different values of q when t is fixed at \(t=10\) and \(\alpha =\beta =0.50\) for Example 1

Fig. 7
figure 7

The graphs of the approximate wave solutions u(xt), v(xt) for Example 2 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=0.25\)

Fig. 8
figure 8

The graphs of the approximate wave solutions u(xt), v(xt) for Example 2 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=0.50\)

Fig. 9
figure 9

The graphs of the approximate wave solutions u(xt), v(xt) for Example 2 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=0.75\)

Fig. 10
figure 10

The graphs of the approximate wave solutions u(xt), v(xt) for Example 2 when \(\alpha =\beta =0.50, \gamma =1\) and \(q=1\)

Fig. 11
figure 11

The graphs of the exact wave solutions u(xt), v(xt) for Example 2 when \(\alpha =\beta =0.50, \gamma =1\)

Fig. 12
figure 12

The graphs of the approximate wave solutions u(xt), v(xt) vs. x among different values of q when t is fixed at \(t=10\) and \(\alpha =\beta =0.50, \gamma =1\) for Example 2

6 Conclusion

In this paper, a new approximate-analytical approach for solving systems of FNPDEs was introduced. The important benefit of this new approximate-analytical approach is that the approximate solutions for fully general systems of FNPDEs can be easily solved without any assumption, and that makes the proposed method much better than the other existing analytical methods. However, we applied the proposed method to obtain the solitary wave solutions and traveling wave solutions for systems of time-fractional dispersive wave equations and the system of the time-fractional long water wave equation. Moreover, the solutions obtained by the proposed method were in excellent agreement with those obtained via previous works and also it was in very good numerical conformity with the exact solutions to confirm the effectiveness and accuracy of the proposed method. The numerical solutions are also carried out in forms of tables and graphs. We used Mathematica and Maple software to obtain the analytical and numerical solutions as well as plotting the graphs.