A Novel Analytical LRPSM for Solving Nonlinear Systems of FPDEs

Abstract

This article employs the Laplace residual power series approach to study nonlinear systems of time-fractional partial differential equations with time-fractional Caputo derivative. The proposed technique is based on a new fractional expansion of the Maclurian series, which provides a rapid convergence series solution where the coefficients of the proposed fractional expansion are computed with the limit concept. The nonlinear systems studied in this work are the Broer-Kaup system, the Burgers’ system of two variables, and the Burgers’ system of three variables, which are used in modeling various nonlinear physical applications such as shock waves, processes of the wave, transportation of vorticity, dispersion in porous media, and hydrodynamic turbulence. The results obtained are reliable, efficient, and accurate with minimal computations. The proposed technique is analyzed by applying it to three attractive problems where the approximate analytical solutions are formulated in rapid convergent fractional Maclurian formulas. The results are studied numerically and graphically to show the performance and validity of the technique, as well as the fractional order impact on the behavior of the solutions. Moreover, numerical comparisons are made with other well-known methods, proving that the results obtained in the proposed technique are much better and the most accurate. Finally, the obtained outcomes and simulation data show that the present method provides a sound methodology and suitable tool for solving such nonlinear systems of time-fractional partial differential equations.

1. Introduction

Fractional-order systems have acquired a lot of attention and interest in various engineering and scientific fields as popular mathematical models used to describe real-world physical phenomena [1,2,3,4,5]. Fractional calculus provides a valuable instrument for showing the development of complicated dynamical systems with long-term memory impacts. In contrast to ordinary derivatives, defining fractional order derivatives of a specific function necessitates the existence of its complete history. Such a non-local feature, i.e., the memory consequence, has made it much more practical to explain various real-world physical systems using fractional differential equations. Investigating dynamics, including complexity, chaos, stability, bifurcation, and synchronization of these fractional order systems, has recently become an interesting research field in nonlinear sciences [6,7,8,9,10,11,12,13]. In order to study the real-world physical systems’ dynamic behavior, it is essential to determine how these solution trajectories can change over slight perturbations. Therefore, performing and developing various numerical techniques to analyze and simulate the systems’ nonlinear dynamics is important. Considering fractional derivatives, analytic-numeric approaches to fractional calculus frequently depend on versions of the Riemann-Liouville, Caputo, Grunwald-Letnikov, Riesz, or other approaches, which were discussed in previous studies during the past few years [14,15,16]. This study, however, will use Caputo’s approach of fractional differentiation, benefiting from Caputo’s approach that initial conditions of the fractional partial differential equations, i.e., (FPDEs) with Caputo’s derivatives take the same conventional form as in integer order.

Differential equations (DEs) can be used for modeling many chemical, biological, and physical phenomena. Because FPDEs have a significant impact on many applied disciplines, particularly nonlinear ones such as fluid flow, biological diffusion of populations, dynamical systems, control theory, electromagnetic waves, etc., there has been a growing interest in them in recent years [17,18,19,20,21]. Most scientific phenomena in various disciplines such as physics, biological systems, and engineering are nonlinear problems; therefore, it might be challenging to find their exact solutions, e.g., physical problems are typically modeled by utilizing higher nonlinear FPDEs, thereby finding exact solutions for these problems is quite challenging. Thus, numerical as well as approximate methods must be employed. Numerous useful techniques were used for solving linear and nonlinear FPDEs, including the variational iteration technique, the Adomian decomposition technique, the homotopy analysis technique, the homotopy perturbation technique, and the fractional residual power series technique [22,23,24,25,26,27,28].

The fractional power series method (FPSM) has been employed to solve several classes of differential and integral equations of the fractional order if the solution of the equation can be extended into a fractional power series [29]. Moreover, FPSM is a fast and easy method utilized to determine the fractional power series solution coefficients because if we compare the computational effort required to compute the solutions of the FPDEs in FPSM with other methods, it becomes clear that it is much less. Moreover, the results are much better, as the speed of implementation on mathematical packages helps to obtain the results in less time and with more accuracy, especially in non-linear problems [30,31,32]. Recently, the FPSM has received the attention of many researchers, whereby various fractional integral and differential equations were investigated successfully by using FPSM, involving fractional Fokker–Planck equations [33], Sawada–Kotera–Ito, Lax, and Kaup–Kupershmidt equations [34], fractional Fredholm integrodifferential equation of order 2β arising in natural sciences [35]. The Laplace transform (LT) technique represents a simple technique for solving several kinds of linear differential integral and integrodifferential equations, as well as a specific class of linear FPDEs [5]. Solving linear DEs by LT technique involves three steps. Transforming the main DEs into the Laplace space represents the first step of this process. Solve the new equation algebraically in the Laplace space in the second step. The last step involves transforming back the obtained solution in the previous step into the initial space, which solves the problem at hand [36].

Overall, there are no semi-approximate or conventional analytical methods that can produce accurate closed-form or approximate solutions for nonlinear FPDE systems. Accordingly, there is a pressing need for efficient numerical methods so that accurate approximate solution can be found for these models for extended periods. Motivated by the above-mentioned discussion, designing an innovative iterative algorithm to produce analytical solutions to the nonlinear FPDE systems is the main aim of our study. The motivation of this study is to present an analytical method called LFPSM to solve a nonlinear system of FPDEs. To specify the efficacy and accuracy of this method, we apply it to solve three nonlinear systems of FPDEs and compare the results obtained with the exact solutions and solutions obtained by other methods. According to our best knowledge, the proposed method has not been applied to find analytical solutions to Broer–Kaup and Burgers’ systems of fractional orders in the literature, which intensely motivated this work.

This study primarily aims to generate accurate approximate solutions to nonlinear FPDE systems in the Caputo sense, which are subject to proper initial conditions by using an innovative analytical algorithm. This algorithm is called Laplace FPSM, which has been suggested and proved in [37]. It is worth mentioning that this newly introduced method relies on transforming the considered equation into the LT space so that a sequence of Laplace series solutions to the new equation form is established, and then the solution to the considered equation can be established by utilizing the inverse LT. Without perturbation, linearization, or discretization, this innovative method can be applied to generate the FPS expansion solutions for both linear and nonlinear FPDEs [38,39]. Furthermore, this technique, unlike the conventional FPSM, does not necessitate matching the corresponding coefficients terms nor the utilization of a relation of recursion. The technique offered is based on the limit concept for finding the variable coefficients. Unlike FPSM, which needs numerous times to compute different fractional derivatives in the steps of the solution, only a few computations are needed to determine the coefficients specified. Therefore, this proposed method has the capability of yielding closed-form solutions, in addition to accurate approximate solutions, by involving a fast convergence series.

The rest of the article is organized as follows. A review of some necessary definitions, properties, and theorems concerning fractional calculus, Laplace transform, and Laplace fractional expansion is presented in Section 2. The methodology for solving a system of nonlinear time-FPDEs by Laplace FPSM is deeply investigated in Section 3. In Section 4, the Broer-Kaup (BK) system of nonlinear time FPDEs, and two Burgers’ systems of nonlinear time FPDEs are solved to show that our approach is accurate and applicable. The results are debated graphically and numerically in Section 5. Finally, Section 6 is lifted for the conclusions.

2. Preliminary Concepts

This section is devoted to overviewing the essential definitions and theorems of fractional differentiation, in addition, to giving a brief for some preliminary definitions and necessary theorems regarding LT, which will be used in sections three and four.

Definition 1.

For$n\in \mathbb{N}$, and$𝒶\in {\mathbb{R}}^{+}$the time-fractional derivative in the Caputo sense for the real-valued function$\mathcal{U}(𝓍,𝓉)$is defined as: [3]

$${\mathfrak{D}}_{𝓉}^{𝒶}\mathcal{U}\left(𝓍,𝓉\right)=\left\{\begin{array}{ll}{\mathcal{I}}_{𝓉}^{n-\alpha }\left(D_t^n\mathcal{U}\left(𝓍,𝓉\right)\right),\hfill & 0<n-1<𝒶\le n, \hfill \\ D_{𝓉}^n\mathcal{U}\left(𝓍,𝓉\right),\hfill & 𝒶=n,\hfill \end{array}\right.$$

where$D_t^n=\frac{{\partial }^n}{\partial t^n},$and${\mathcal{I}}_{\mathit{t}}$is the R-L fractional integral operator and which is given by:

$${\mathcal{I}}_{𝓉}^{𝒶}\mathcal{U}\left(𝓍,𝓉\right)=\left\{\begin{array}{ll}\frac{1}{\mathsf{\Gamma }\left(𝒶\right)}\underset{0}{\overset{t}{{\displaystyle \int }}}\frac{\mathcal{U}\left(x,\eta \right)}{{\left(𝓉-\eta \right)}^{1-𝒶}}d\eta , \hfill & 0\le \eta <𝓍,𝓉>0,\hfill \\ \mathcal{U}\left(𝓍,𝓉\right)\hfill & 𝒶=0.\hfill \end{array}\right.$$

Consequently, for$n-1<𝒶\le n, \beta >-1$and$𝓉\ge 0$, the operators${\mathfrak{D}}_{𝓉}^{𝒶}$and${\mathcal{I}}_{𝓉}^{\alpha }$satisfy the following properties:

• ${\mathfrak{D}}_{𝓉}^{𝒶}c=0$, $c\in ℝ$.
• ${\mathfrak{D}}_{𝓉}^{𝒶}{𝓉}^{\beta }=\frac{\mathsf{\Gamma }\left(\beta +1\right)}{\mathsf{\Gamma }\left(\beta +1-𝒶\right)}{𝓉}^{\beta -𝒶}.$
• ${\mathfrak{D}}_{𝓉}^{𝒶}{\mathcal{I}}_{𝓉}^{𝒶}\mathcal{U}\left(𝓍,𝓉\right)=\mathcal{U}\left(𝓍,𝓉\right)$.
• ${\mathcal{I}}_{𝓉}^{𝒶}{\mathfrak{D}}_{𝓉}^{𝒶}\mathcal{U}\left(𝓍,𝓉\right)=\mathcal{U}\left(𝓍,𝓉\right)-\sum _{j=0}^{n-1}{D}_{𝓉}^j\left(𝓍,0^{+}\right)\frac{{𝓉}^j}{j!}$, for$\mathcal{U}\in C^n\left[a,b\right]$, n$-1<𝒶\le n, n\in \mathbb{N} and a,b\in ℝ$.

Definition 2.

The Laplace transformation (LT) of the piecewise continuous function $\mathcal{U}\left(𝓍,𝓉\right)$on$I\times \left[0\right.,\left.\infty \right)$and of exponential order$\delta$is given by: [38]

$$\mathfrak{U}\left(𝓍,\mathfrak{s}\right)=\mathcal{L}\left[\mathcal{U}\left(𝓍,𝓉\right)\right]≔{{\displaystyle \int }}_0^{\infty }{e}^{-\mathfrak{s}𝓉}\mathcal{U}\left(𝓍,𝓉\right)d𝓉, \mathfrak{s}>\delta ,$$

and the inverse LT of the transform function$\mathfrak{U}\left(𝓍,\mathfrak{s}\right)$is given by:

$$\mathcal{U}\left(𝓍,𝓉\right)={\mathcal{L}}^{-1}\left[\mathfrak{U}\left(𝓍,\mathfrak{s}\right)\right]={{\displaystyle \int }}_{c-i\infty }^{c+i\infty }{e}^{\mathfrak{s}𝓉}\mathfrak{U}\left(𝓍,\mathfrak{s}\right)d\mathfrak{s},c=Re\left(\mathfrak{s}\right)>{\delta }_0,$$

where${\delta }_0$lies in the right half plane of the absolute convergence of the Laplace integral.

Lemma 1.

Let$\mathcal{U}\left(𝓍,𝓉\right)$and$\mathcal{V}\left(𝓍,𝓉\right)$be two piecewise continuous functions defined on$I\times \left[0\right.,\left.\infty \right)$and of exponential order${\delta }_1$and${\delta }_2$, respectively, where${\delta }_1<{\delta }_2$. Suppose that$\mathfrak{U}\left(𝓍,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{U}\left(𝓍,𝓉\right)\right], \bigvee \left(𝓍,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{V}\left(𝓍,𝓉\right)\right]$and$\left\{a,b\right\}\in ℝ$. Then, [38]

• $\mathcal{L}\left[a\mathcal{U}\left(𝓍,𝓉\right)+b\mathcal{V}\left(𝓍,𝓉\right)\right]=a\mathfrak{U}\left(𝓍,\mathfrak{s}\right)+b\bigvee \left(𝓍,\mathfrak{s}\right),𝓍\in I, \mathfrak{s}>{\delta }_1$.
• ${\mathcal{L}}^{-1}\left[a\mathfrak{U}\left(𝓍,\mathfrak{s}\right)+bx\left(\bigvee ,\mathfrak{s}\right)\right]=a\mathcal{U}\left(𝓍,𝓉\right)+b\mathcal{V}\left(𝓍,𝓉\right),x\in I, 𝓉\ge 0$.
• $\mathcal{L}\left[e^{a𝓉}\mathcal{U}\left(𝓍,𝓉\right)\right]=\mathfrak{U}\left(𝓍,\mathfrak{s}-a\right),x\in I, \mathfrak{s}>a+{\delta }_1$.
• $\underset{\mathfrak{s}\to \infty }{\lim }\mathfrak{s}\mathfrak{U}\left(𝓍,\mathfrak{s}\right)=\mathcal{U}\left(𝓍,0\right), 𝓍\in I$.

Lemma 2.

Let$\mathcal{U}\left(𝓍,𝓉\right)$be a piecewise continuous function defined on$I\times \left[0\right.,\left.\infty \right)$and of exponential order$\delta$, and$\mathfrak{U}\left(𝓍,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{U}\left(𝓍,𝓉\right)\right]$. Then, [31]

• $\mathcal{L}\left[\mathcal{I}\mathcal{U}\left(𝓍,𝓉\right)\right]={\mathfrak{s}}^{𝒶-1}\mathfrak{U}\left(x𝓍,\mathfrak{s}\right), 𝒶>0$.
• $\mathcal{L}\left[{\mathfrak{D}}_{𝓉}^{𝒶}\mathcal{U}\left(𝓍,𝓉\right) \right]={\mathfrak{s}}^{𝒶}\mathfrak{U}\left(𝓍,\mathfrak{s}\right)-{\displaystyle \sum }_{k=0}^{n-1}{\mathfrak{s}}^{𝒶-k-1}{D}_t^k\mathcal{U}\left(𝓍,0\right),n-1< 𝒶\le n$.
• $\mathcal{L}\left[{\mathfrak{D}}_{𝓉}^{j𝒶}\mathcal{U}\left(𝓍,𝓉\right) \right]={\mathfrak{s}}^{j𝒶}\mathfrak{U}\left(𝓍,\mathfrak{s}\right)-{\displaystyle \sum }_{k=0}^{j-1}{\mathfrak{s}}^{\left(j-k\right)𝒶-1}{D}_{𝓉}^{k𝒶}\mathcal{U}\left(𝓍,0\right),0<𝒶\le 1,$where${\mathfrak{D}}_{𝓉}^{j𝒶}={\mathfrak{D}}_{𝓉}^{𝒶}·{\mathfrak{D}}_{𝓉}^{𝒶}···{\mathfrak{D}}_{𝓉}^{𝒶}$(j-times).

Proof. 

The proof is in [38]. □

Theorem 1.

Let $\mathcal{U}\left(𝓍,𝓉\right)$be a piecewise continuous function defined on$I\times \left[0\right.,\left.\infty \right)$and of exponential order$\delta$. Suppose that the function$\mathfrak{U}\left(𝓍,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{U}\left(𝓍,𝓉\right)\right]$has the following fractional expansion (FE): [38]

$$\mathfrak{U}\left(𝓍,\mathfrak{s}\right)={\displaystyle \sum }_{n=0}^{\infty }\frac{{𝒽}_n\left(𝓍\right)}{{\mathfrak{s}}^{n\alpha +1}}, x\in I, \mathfrak{s}>\delta , 0<𝒶\le 1.$$

Then,${𝒽}_n\left(𝓍\right)={\mathfrak{D}}_{𝓉}^{n𝒶}\mathcal{U}\left(𝓍,0\right)$.

Proof. 

The proof is in [38]. □

Remark 1.

The inverse LT${\mathcal{L}}^{-1}\left[\mathfrak{U}\left(𝓍,\mathfrak{s}\right)\right]=\mathcal{U}\left(𝓍,𝓉\right),$in Theorem 1 is in the following expansion series (FSE) form:

$$\mathcal{U}\left(𝓍,𝓉\right)={\displaystyle \sum }_{n=0}^{\infty }{\mathfrak{D}}_{𝓉}^{n𝒶}\mathcal{U}\left(𝓍,0\right)\frac{{𝓉}^{n𝒶}}{\mathsf{\Gamma }\left(n𝒶+1\right)}, 0< \alpha \le 1,𝓉>0.$$

Theorem 2.

Let$\mathcal{U}\left(𝓍,𝓉\right)$be an exponential function of order$\delta$defined on$I\times \left[0\right.,\left.\infty \right)$, and let$\mathfrak{U}\left(𝓍,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{U}\left(𝓍,𝓉\right)\right]$can be represented as the FE in Theorem 1. If$\left|\mathfrak{s}\mathcal{L}\left[{\mathfrak{D}}_{𝓉}^{\left(n+1\right)𝒶}\mathcal{U}\left(𝓍,𝓉\right) \right]\right|\le M\left(𝓍\right)$, on$I\times \left(\delta \right.,\left.\gamma \right]$where$0< 𝒶\le 1$, then the reminder$R_n\left(𝓍,\mathfrak{s}\right)$of the FE in Theorem 1 satisfies the following inequality: [38]

$$\left|R_n\left(𝓍,\mathfrak{s}\right)\right|\le \frac{M\left(𝓍\right)}{{\mathfrak{s}}^{1+\left(n+1\right)𝒶}}, x\in I, \delta <\mathfrak{s}\le \gamma$$

Proof. 

The proof is in [38]. □

Theorem 3.

If$𝒶\in \left.0,1\right]$,$\parallel {\mathcal{U}}_{k+1}\left(𝓍,𝓉\right)\parallel \le 𝒶\parallel {\mathcal{U}}_k\left(𝓍,𝓉\right)\parallel$gives$\forall k\in N$and$0<t<{\rm T}<1$, then the series of numerical solutions converges to the exact solution [39].

Proof. 

We notice that $\forall 0<𝓉<\mathsf{{\rm T}}<1$,

$$\begin{array}{c}\parallel \mathcal{U}\left(𝓍,𝓉\right)-{\mathcal{U}}_k\left(𝓍,𝓉\right)\parallel =\parallel {\displaystyle {\displaystyle \sum }_{m=k+1}^{\infty }}{\mathcal{U}}_m\left(𝓍,𝓉\right)\parallel \le {\displaystyle {\displaystyle \sum }_{m=k+1}^{\infty }}\parallel {\mathcal{U}}_m(𝓍,𝓉\parallel )\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le \parallel ℊ\left(\eta \right)\parallel \parallel {\displaystyle {\displaystyle \sum }_{m=k+1}^{\infty }}{C}^m\parallel =\frac{{𝒶}^{k+1}}{1-𝒶}\parallel ℊ\left(\eta \right)\parallel \to 0 \mathrm{as} k\to \infty .\square \end{array}$$

3. The Methodology of Laplace RPSM

In this part, we present the fundamental idea of the Laplace RPSM for solving the system of time FPDEs with initial conditions. Our strategy for using the proposed scheme is to rely on coupling the Laplace transform and the RPS approach. More precisely, consider the system of FPDEs with the initial conditions of the form:

$$\left\{\begin{array}{c}{\mathfrak{D}}_t^{𝒶}\mathit{U}\left(\mathbf{\eta },𝓉\right)={\mathcal{A}}_1\left[\mathit{U}\left(\mathbf{\eta },𝓉\right)\right]+{\mathcal{A}}_2\left[\mathit{U}\left(\mathbf{\eta },𝓉\right)\right], 0<𝒶\le 1,\\ \hspace{1em}\mathit{U}\left(\mathbf{\eta },0\right)={\mathcal{P}}_j\left(\mathbf{\eta }\right), j=1,2,\dots ,n ,\end{array}\right.$$

(1)

where ${\mathcal{A}}_1$, ${\mathcal{A}}_2$ are two linear or nonlinear operators such that $\mathit{U}\left(\mathbf{\eta },𝓉\right)=\left({\mathcal{U}}_1\left(\mathbf{\eta },𝓉\right),{\mathcal{U}}_2\left(\mathbf{\eta },𝓉\right),\dots ,{\mathcal{U}}_n\left(\mathbf{\eta },𝓉\right)\right),$ is the unknown vector function to be determined, and $\mathbf{\eta }= \left({\eta }_1, {\eta }_2,\dots ,{\eta }_m\right)\in {ℝ}^m, n,m\in ℝ$. Here, ${\mathfrak{D}}_t^{𝒶}$ refers to the time-fractional derivative of order $𝒶\in \left(0,1\right]$, in the Caputo meaning.

To build the approximate solution of (1) by using the Laplace RPSM, one can accomplish the following procedure:

Step 1: Taking the LT on the two sides of (1) and employing the initial data of (1), as well as relying on Lemma 2, part (2), we get:

$$\begin{gathered} \mathit{u}\left(\mathbf{\eta },\mathfrak{s}\right)=\frac{{\mathcal{P}}_j\left(\mathbf{\eta }\right)}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\mathcal{L}\left\{{\mathcal{A}}_1\left[\mathit{U}\left(\mathbf{\eta },t\right)\right]\}+\mathcal{L}\{{\mathcal{A}}_2\left[\mathit{U}\left(\mathbf{\eta },𝓉\right)\right]\right\}\right), \\ \mathrm{where} \mathit{u}\left(\mathbf{\eta },\mathfrak{s}\right)= \mathcal{L}\left[\mathit{U}\left(\mathbf{\eta },𝓉\right)\right]\left(\mathfrak{s}\right), \mathfrak{s}>\delta . \end{gathered}$$

(2)

Step 2: Based on Theorem 1, we suppose that the approximate solution of the Laplace Equation (2) has the following Laplace fractional expansions:

$$\begin{gathered} {𝓊}_1\left(\mathbf{\eta },\mathfrak{s}\right)=\frac{{\mathcal{P}}_1\left(\mathbf{\eta }\right)}{\mathfrak{s}}+{\displaystyle \sum }_{n=1}^{\infty }\frac{{𝒽}_n\left(\mathbf{\eta }\right)}{{\mathfrak{s}}^{n𝒶+1}}, \mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0, \\ {𝓊}_2\left(\mathbf{\eta },\mathfrak{s}\right)=\frac{{\mathcal{P}}_2\left(\mathbf{\eta }\right)}{\mathfrak{s}}+{\displaystyle \sum }_{n=1}^{\infty }\frac{{𝒽}_n\left(\mathbf{\eta }\right)}{{\mathfrak{s}}^{n\mathbf{𝒶}+1}}, \mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0, \\ ⋮ \\ {𝓊}_n\left(\mathbf{\eta },\mathfrak{s}\right)=\frac{{\mathcal{P}}_n\left(\mathbf{\eta }\right)}{\mathfrak{s}}+{\displaystyle \sum }_{n=1}^{\infty }\frac{{𝒽}_n\left(\mathbf{\eta }\right)}{{\mathfrak{s}}^{n\mathbf{𝒶}+1}}, \mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0, \end{gathered}$$

(3)

and the $k-th$ Laplace series solutions take the following form:

$$\begin{gathered} u_{1,k}\left(\mathbf{\eta },\mathfrak{s}\right)=\frac{{\mathcal{P}}_1\left(\mathbf{\eta }\right)}{\mathfrak{s}}+{\displaystyle \sum }_{n=1}^k\frac{{𝒽}_n\left(\mathbf{\eta }\right)}{{\mathfrak{s}}^{n𝒶+1}}, \mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0, \\ {u}_{2,k}\left(\mathbf{\eta },\mathfrak{s}\right)=\frac{{\mathcal{P}}_2\left(\mathbf{\eta }\right)}{\mathfrak{s}}+{\displaystyle \sum }_{n=1}^k\frac{{𝒽}_n\left(\mathbf{\eta }\right)}{{\mathfrak{s}}^{n𝒶+1}}, \mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0, \\ ⋮ \\ {u}_{n,k}\left(\mathbf{\eta },\mathfrak{s}\right)=\frac{{\mathcal{P}}_n\left(\mathbf{\eta }\right)}{\mathfrak{s}}+{\displaystyle \sum }_{n=1}^k\frac{{𝒽}_n\left(\mathbf{\eta }\right)}{{\mathfrak{s}}^{n𝒶+1}}, \mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0. \end{gathered}$$

(4)

Step 3: Define the $k-th$ Laplace fractional residual function of (2) as:

$$\mathcal{L}\left(Re{s}_{u_k}\left(\mathbf{\eta },\mathfrak{s}\right)\right)=\frac{{\mathcal{P}}_j\left(\mathbf{\eta }\right)}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\mathcal{L}\left\{{\mathcal{A}}_1\left[\mathit{U}\left(\mathbf{\eta },𝓉\right)\right]\}+\mathcal{L}\{{\mathcal{A}}_2\left[\mathit{U}\left(\mathbf{\eta },𝓉\right)\right]\right\}\right),$$

(5)

and the Laplace fractional residual function of (2) can be defined as:

$$\begin{gathered} \underset{k\to \infty }{\lim }\mathcal{L}\left(Re{s}_{{\mathbf{𝓾}}_k}\left(\mathbf{\eta },\mathfrak{s}\right)\right)=\mathcal{L}\left(Re{s}_{\mathbf{𝓾}}\left(\mathbf{\eta },\mathfrak{s}\right)\right) \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{{\mathcal{P}}_j\left(\mathbf{\eta }\right)}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\mathcal{L}\left\{{\mathcal{A}}_1\left[\mathit{U}\left(\mathbf{\eta },t\right)\right]\}+\mathcal{L}\{{\mathcal{A}}_2\left[\mathit{U}\left(\mathbf{\eta },t\right)\right]\right\}\right). \end{gathered}$$

(6)

As in [37,38,39], some of the beneficial facts of Laplace residual function, which are fundamental in constructing the approximate solution, are listed as follows:

• $\underset{k\to \infty }{lim}\mathcal{L}\left(Re{s}_{{\mathbf{𝓾}}_k}\left(\mathbf{\eta },\mathfrak{s}\right)\right)=\mathcal{L}\left(Re{s}_{\mathbf{𝓾}}\left(\mathbf{\eta },\mathfrak{s}\right)\right),$ for $\mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0$.
• $\mathcal{L}\left(Re{s}_{\mathbf{𝓾}}\left(\mathbf{\eta },\mathfrak{s}\right)\right)=0$, for $\mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0$.
• $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{k𝒶+1}\mathcal{L}\left(Re{s}_{{\mathbf{𝓾}}_k}\left(\mathbf{\eta },\mathfrak{s}\right)\right)=0$, for $\mathbf{\eta }\in \overline{I}, \mathfrak{s}>\delta \ge 0$, and $k=1,2,3,\dots$

Step 4: The $k-th$ Laplace fractional residual function of (5) is substituted by the $k-th$ Laplace series solution (4).

Step 5: By solving the system $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{k𝒶+1}\mathcal{L}\left(Re{s}_{u_k}\left(\mathbf{\eta },\mathfrak{s}\right)\right)=0$, the unknown coefficients ${𝒽}_k\left(\mathbf{\eta }\right)$, for $k=1,2,3,\dots$, easily could be founded. Then, we accumulate the received variable coefficients in terms of the Laplace fractional expansion series (4) ${\mathit{u}}_{j,k}\left(\mathbf{\eta },\mathfrak{s}\right)$.

Step 6: The approximate solution ${\mathit{U}}_{\mathit{j},\mathit{k}}\left(\mathbf{\eta },t\right)$, of the main Equation (1), can be attained by applying the inverse Laplace transform operator on both sides of the obtained Laplace series solution.

4. Numerical Examples

In this section, we show that the Laplace RPSM is superior, efficient, and applicable, which is achieved by testing three nonlinear time-FPDEs systems. It should be noted here that all numerical and symbolic calculations are made using the Mathematica 12 software package.

Example 1.

Consider the following Broer-Kaup system of nonlinear time-FPDEs:

$$\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{{\partial }^{𝒶}\mathcal{U}}{\partial {𝓉}^{𝒶}}+\mathcal{U}\frac{\partial \mathcal{U}}{\partial x}+\frac{\partial \mathcal{V}}{\partial x}=0,\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{{\partial }^{𝒶}\mathcal{V}}{\partial {𝓉}^{𝒶}}+\frac{\partial \mathcal{U}}{\partial x}+\frac{\partial \left(\mathcal{U}\mathcal{V}\right)}{\partial x}+\frac{{\partial }^3\mathcal{U}}{\partial x^3}=0,\\ \mathit{subject}\mathit{to}\mathit{ICs}& \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathcal{U}\left(x,0\right)=1+2\mathit{tanh}\left(x\right), \mathrm{and} \mathcal{V}\left(x,0\right)=1-2{\mathit{tanh}}^2\left(x\right), \end{array}$$

(7)

where $𝒶\in \left(0,1\right]$and$\left(x,t\right)\in ℝ\times \left[0,1\right]$. The exact solutions when$𝒶=1$, are$\left(\mathcal{U}\left(x,t\right),\mathcal{V}\left(x,t\right)\right)=\left(1-2\mathit{tanh}\left(t-x\right),1-2{\mathit{tanh}}^2\left(x-t\right)\right)$.

By applying the LT operator on (7) and using the second part of Lemma 2 and the ICs of (7), the Laplace fractional equations are:

$$\begin{array}{l}\mathfrak{U}\left(x,\mathfrak{s}\right)=\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\mathrm{\bigvee }\right\}\right\},\\ \hspace{1em}\mathrm{\bigvee }\left(x,\mathfrak{s}\right)=\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\right\}-\mathcal{L}\frac{\partial }{\partial x}\left\{{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}{\mathcal{L}}^{-1}\left\{v\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^3}{\partial x^3}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\right\},\end{array}$$

(8)

where $\mathfrak{U}\left(x,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{U}\left(x,t\right)\right]$ and $\mathrm{\bigvee }\left(x,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{V}\left(x,t\right)\right]$.

According to the last discussion of the proposed method, the $k-th$ Laplace series solutions, ${\mathfrak{U}}_k\left(x,\mathfrak{s}\right)$ and ${\bigvee }_k\left(x,\mathfrak{s}\right)$ for (8) are expressed as:

$$\begin{array}{l}{\mathfrak{U}}_k\left(x,\mathfrak{s}\right)=\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{𝒽}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}},\\ \hspace{1em}{\bigvee }_k\left(x,\mathfrak{s}\right)=\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{ℊ}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}}.\end{array}$$

(9)

Hence, the $k-th$ Laplace fractional residual functions of (8) is defined as:

$$\begin{gathered} \mathcal{L}\left(Re{s}_{{\mathfrak{U}}_k}\left(x,\mathfrak{s}\right)\right)={\displaystyle \sum }_{n=1}^k\frac{{𝒽}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\right\}, \\ \mathcal{L}\left(Re{s}_{{\bigvee }_k}\left(x,\mathfrak{s}\right)\right)={\displaystyle \sum }_{n=1}^k\frac{{ℊ}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\right\}+\mathcal{L}\frac{\partial }{\partial x}\left\{{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^3}{\partial x^3}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\right\}. \end{gathered}$$

(10)

The $1-st$ Laplace fractional residual functions can be carried out by letting $k=1$, in (10):

$$\begin{array}{cc}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_1}\left(x,\mathfrak{s}\right)\right)=\hfill & \frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & =\frac{1}{{\mathfrak{s}}^{𝒶+1}}\left({𝒽}_1\left(x\right)+2{\sec \mathrm{h}}^2\left(x\right)\right)\hfill \\ & +\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left(2{𝒽}_1\left(x\right){\sec \mathrm{h}}^2\left(x\right)+{{ℊ}^{\prime }}_1\left(x\right)+{{𝒽}^{\prime }}_1\left(x\right)+2{{𝒽}^{\prime }}_1\left(x\right)\tanh \left(x\right)\right)\hfill \\ & +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left({𝒽}_1\left(x\right){{𝒽}^{\prime }}_1\left(x\right)\right)\frac{\mathsf{\Gamma }\left(2𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(𝒶+1\right)\right)}^2},\hfill \\ \mathcal{L}\left(Re{s}_{{\bigvee }_1}\left(x,\mathfrak{s}\right)\right)=\hfill & \frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\frac{\partial }{\partial x}\left\{{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}{\mathcal{L}}^{-1}\left\{\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^3}{\partial x^3}{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & =\frac{1}{{\mathfrak{s}}^{𝒶+1}}\left({ℊ}_1\left(x\right)-4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)\right)\hfill \\ & +\frac{1}{{\mathfrak{s}}^{2𝒶+1}}(2{ℊ}_1\left(x\right){\sec \mathrm{h}}^2\left(x\right)-4{𝒽}_1\left(x\right)\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)+{{ℊ}^{\prime }}_1\left(x\right)+2{{ℊ}^{\prime }}_1\left(x\right)\tanh \left(x\right)\hfill \\ & +2{{𝒽}^{\prime }}_1\left(x\right){\sec \mathrm{h}}^2\left(x\right)+{𝒽}_1{}^{\left(3\right)}\left(x\right))+\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left({𝒽}_1\left(x\right){{ℊ}^{\prime }}_1\left(x\right)+{ℊ}_1\left(x\right){{𝒽}^{\prime }}_1\left(x\right)\right)\frac{\mathsf{\Gamma }\left(2𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(𝒶+1\right)\right)}^2}.\hfill \end{array}$$

(11)

To find the $1-st$ Laplace series solution of (8), we simply take the next process $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{x+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_1}\left(x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\bigvee }_1}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right)$, which yields that ${𝒽}_1\left(x\right)=-2{\sec \mathrm{h}}^2\left(x\right)$ and ${ℊ}_1\left(x\right)=4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)$. So, the $1-st$ Laplace series solutions of (8) are:

$$\begin{array}{l}{\mathfrak{U}}_1\left(x,\mathfrak{s}\right)=\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}},\\ {\bigvee }_1\left(x,\mathfrak{s}\right)=\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}.\end{array}$$

(12)

For $k=2$, in (10) the $2-nd$ Laplace residual functions can be written as:

$$\begin{array}{cc}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_2}\left(x,\mathfrak{s}\right)\right)=\hfill & -\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\right\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\hfill \\ & +\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ & =\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left({𝒽}_2\left(x\right)+4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)\right)\hfill \\ & +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}(2{𝒽}_2\left(x\right){\sec \mathrm{h}}^2\left(x\right)+{{ℊ}^{\prime }}_2\left(x\right)+{{𝒽}^{\prime }}_2\left(x\right)\hfill \\ & +2{{𝒽}^{\prime }}_2\left(x\right)\tanh \left(x\right)-8\tanh \left(x\right){\sec \mathrm{h}}^4\left(x\right)\frac{\mathsf{\Gamma }\left(2𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(𝒶+1\right)\right)}^2})\hfill \\ & +\frac{1}{{\mathfrak{s}}^{4𝒶+1}}\left(4{𝒽}_2\left(x\right)\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)-2{{𝒽}^{\prime }}_2\left(x\right){\sec \mathrm{h}}^2\left(x\right)\right)\frac{\mathsf{\Gamma }\left(3𝒶+1\right)}{\mathsf{\Gamma }\left(2𝒶+1\right)\mathsf{\Gamma }\left(𝒶+1\right)}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{5𝒶+1}}\left({𝒽}_2\left(x\right){{𝒽}^{\prime }}_2\left(x\right)\right)\frac{\mathsf{\Gamma }\left(4𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(2𝒶+1\right)\right)}^2},\hfill \\ \mathcal{L}\left(Re{s}_{{\bigvee }_2}\left(x,\mathfrak{s}\right)\right)=\hfill & \frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\frac{\partial }{\partial x}\left\{{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}{\mathcal{L}}^{-1}\right\{\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}\hfill \\ & +\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^3}{\partial x^3}{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ & =\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left({ℊ}_2\left(x\right)-4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)\right)\hfill \\ & +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}({𝒽}_2{}^{\left(3\right)}\left(x\right)-4{𝒽}_2\left(x\right)\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)+{{ℊ}^{\prime }}_2\left(x\right)+2{{ℊ}^{\prime }}_2\left(x\right)\tanh \left(x\right)\hfill \\ & +2{{𝒽}^{\prime }}_2\left(x\right){\sec \mathrm{h}}^2\left(x\right)+2{ℊ}_2\left(x\right){\sec \mathrm{h}}^2\left(x\right)\hfill \\ & +\left(\left(16\cosh \left(2x\right){\sec \mathrm{h}}^6\left(x\right)-24{\sec \mathrm{h}}^6\left(x\right)\right)\frac{\mathsf{\Gamma }\left(2𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(𝒶+1\right)\right)}^2})\right)\hfill \\ & +\frac{1}{{\mathfrak{s}}^{4𝒶+1}}(4{ℊ}_2\left(x\right)\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)-8{𝒽}_2\left(x\right){\sec \mathrm{h}}^2\left(x\right)+12{𝒽}_2\left(x\right){\sec \mathrm{h}}^4\left(x\right)\hfill \\ & -2{{ℊ}^{\prime }}_2\left(x\right){\sec \mathrm{h}}^2\left(x\right)+4{{𝒽}^{\prime }}_2\left(x\right)\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right))\frac{\mathsf{\Gamma }\left(3𝒶+1\right)}{\mathsf{\Gamma }\left(2𝒶+1\right)\mathsf{\Gamma }\left(𝒶+1\right)}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{5𝒶+1}}\left({𝒽}_2\left(x\right)ℊ{\prime }_2\left(x\right)+{ℊ}_2\left(x\right)𝒽{\prime }_2\left(x\right)\right)\frac{\mathsf{\Gamma }\left(4𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(2𝒶+1\right)\right)}^2}.\hfill \end{array}$$

(13)

To find the $2-nd$ Laplace series solution of (8), we simply find out the next process $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{2𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_2}\left(x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{x_2}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right),$ and by solving limits, we get ${𝒽}_2\left(x\right)=-4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)$ and ${ℊ}_2\left(x\right)=4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)$. So, the $2-nd$ Laplace series solution of (8) could be expressed as:

$$\begin{gathered} {\mathfrak{U}}_2\left(x,\mathfrak{s}\right)=\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}, \\ {\bigvee }_2\left(x,\mathfrak{s}\right)=\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)}{{\mathfrak{s}}^{2𝒶+1}}. \end{gathered}$$

(14)

Similarly, for $k=3$, we have:

$$\begin{array}{cc}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_3}\left(x,\mathfrak{s}\right)\right)=\hfill & -\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\right\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}\hfill \\ & -\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\left\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\right\{\frac{1-2\tan {\mathrm{h}}^2\left(x\right)}{\mathfrak{s}}+\frac{4\tan \mathrm{h}\left(x\right)\sec {\mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{4\sec {\mathrm{h}}^4\left(x\right)(\cosh \left(2,x\right)-2)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\left\}\right\}\hfill \\ \mathcal{L}\left(Re{s}_{{\bigvee }_3}\left(x,\mathfrak{s}\right)\right)=\hfill & \frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ & +\mathcal{L}\frac{\partial }{\partial x}\left\{{\mathcal{L}}^{-1}\right\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\hfill \\ & +\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\left\}{\mathcal{L}}^{-1}\left\{\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^3}{\partial x^3}{\mathcal{L}}^{-1}\left\{\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}.\hfill \end{array}$$

(15)

By solving $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{3𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_3}\left(x,\mathfrak{s}\right)\right), \mathcal{L}\left(Re{s}_{{\bigvee }_3}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right)$. It yields that: ${𝒽}_3\left(x\right)=-4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)$ and ${ℊ}_3\left(x\right)=8{\sec \mathrm{h}}^4\left(x\right)\tanh \left(x\right)\left(\cosh \left(2x\right)-5\right)$. So, the $3-rd$ Laplace series solution of (8) could be written as:

$$\begin{gathered} {\mathfrak{U}}_3\left(x,\mathfrak{s}\right)=\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{2{\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}-\frac{4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)}{{\mathfrak{s}}^{3𝒶+1}}, \\ {\bigvee }_3\left(x,\mathfrak{s}\right)=\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+\frac{4\tanh \left(x\right){\sec \mathrm{h}}^2\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{4{\sec \mathrm{h}}^4\left(x\right)\left(\cosh \left(2x\right)-2\right)}{{\mathfrak{s}}^{2𝒶+1}} \\ +\frac{8{\sec \mathrm{h}}^4\left(x\right)\tanh \left(x\right)\left(\cosh \left(2x\right)-5\right)}{{\mathfrak{s}}^{3𝒶+1}}. \end{gathered}$$

(16)

Using Mathematica, we can perform the aforesaid steps for an arbitrary $k$, and using the fact $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{k𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_k}\left(x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\bigvee }_k}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right)$, one can obtain that ${𝒽}_k\left(x\right)={\left(-1\right)}^k\frac{d^{\left(k\right)}}{d{x}^{\left(k\right)}}\left(2\tanh \left(x\right)\right)$ and ${ℊ}_k\left(x\right)={\left(-1\right)}^k\frac{d^{\left(k\right)}}{d{x}^{\left(k\right)}}\left(-2{\tanh }^2\left(x\right)\right)$. Thus, the $k-th$ Laplace series solution of (8) could be reformulated by the following fractional expansions:

$$\begin{array}{cc}{\mathfrak{U}}_k\left(x,\mathfrak{s}\right)=\hfill & (\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}-\frac{\frac{d}{dx}\left(2\tanh \left(x\right)\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\frac{d^{\left(2\right)}}{d{x}^{\left(2\right)}}\left(2\tanh \left(x\right)\right)}{{\mathfrak{s}}^{2𝒶+1}}-\frac{\frac{d^{\left(3\right)}}{d{x}^{\left(3\right)}}\left(2\tanh \left(x\right)\right)}{{\mathfrak{s}}^{3𝒶+1}}+\dots \hfill \\ & +{\left(-1\right)}^k\frac{\frac{d^{\left(k\right)}}{d{x}^{\left(k\right)}}\left(2\tanh \left(x\right)\right)}{{\mathfrak{s}}^{k𝒶+1}})=\frac{1+2\tanh \left(x\right)}{\mathfrak{s}}{\displaystyle {\displaystyle \sum }_{n=1}^k}{\left(-1\right)}^n\frac{\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(2\tanh \left(x\right)\right)}{{\mathfrak{s}}^{n𝒶+1}},\hfill \\ {\bigvee }_k\left(x,\mathfrak{s}\right)=\hfill & (\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}-\frac{\frac{d}{dx}\left(-2{\tanh }^2\left(x\right)\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\frac{d^{\left(2\right)}}{d{x}^{\left(2\right)}}\left(-2{\tanh }^2\left(x\right)\right)}{{\mathfrak{s}}^{2𝒶+1}}-\frac{\frac{d^{\left(3\right)}}{d{x}^{\left(3\right)}}\left(-2{\tanh }^2\left(x\right)\right)}{{\mathfrak{s}}^{3𝒶+1}}+\dots \hfill \\ & +{\left(-1\right)}^k\frac{\frac{d^{\left(k\right)}}{d{x}^{\left(k\right)}}\left(-2{\tanh }^2\left(x\right)\right)}{{\mathfrak{s}}^{k𝒶+1}})=\frac{1-2{\tanh }^2\left(x\right)}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}{\left(-1\right)}^n\frac{\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(-2{\tanh }^2\left(x\right)\right)}{{\mathfrak{s}}^{n𝒶+1}}.\hfill \end{array}$$

(17)

Finally, by applying the inverse Laplace transform for the obtained expansions (17), we conclude that the $k-th$ approximate solution of the time-fractional nonlinear system (7) can be formulated as:

$$\begin{array}{l}{\mathcal{U}}_k\left(x,x\right)=1+2\tanh \left(x\right)+{\displaystyle {\displaystyle \sum }_{n=1}^k}{\left(-1\right)}^n\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(2\tanh \left(x\right)\right)\frac{x^{nx}}{\mathsf{\Gamma }\left(n𝒶+1\right)},\\ \hspace{1em}{\mathcal{V}}_k\left(x,x\right)=1-2{\tanh }^2\left(x\right)+{\displaystyle {\displaystyle \sum }_{n=1}^k}{\left(-1\right)}^n\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(-2{\tanh }^2\left(x\right)\right)\frac{x^{nx}}{\mathsf{\Gamma }\left(n𝒶+1\right)}.\end{array}$$

(18)

When $k\to \infty$ and $𝒶=1$ in (18), we obtain the Maclaurin series expansions of the closed form:

$$\begin{array}{c}\mathcal{U}\left(x,𝓉\right)=1+2\tanh \left(x\right)+{\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{\left(-1\right)}^n\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(2\tanh \left(x\right)\right)\frac{{𝓉}^n}{n!}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=1\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}{\left(-1\right)}^n\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(2\tanh \left(x\right)\right)\frac{{𝓉}^n}{n!}=1+{\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}{\left(-1\right)}^n\left(\frac{d^{\left(n\right)}}{d{𝓉}^{\left(n\right)}}{\left(2\tanh \left(x+𝓉\right)\right)}_{\left|𝓉=0\right.}\right)\frac{{𝓉}^n}{n!}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=1+{\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}\left(\frac{d^{\left(n\right)}}{d{𝓉}^{\left(n\right)}}{\left(2\tanh \left(x-𝓉\right)\right)}_{\left|𝓉=0\right.}\right)\frac{{𝓉}^n}{n!}=1+2\tanh \left(x-𝓉\right)=1-2\tanh \left(𝓉-x\right),\hfill \\ \hspace{1em}\hspace{1em}\mathcal{V}\left(x,x\right)=1-2{\tanh }^2\left(x\right)+{\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{\left(-1\right)}^n\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(-2{\tanh }^2\left(x\right)\right)\frac{{𝓉}^n}{n!}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=1\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}{\left(-1\right)}^n\frac{d^{\left(n\right)}}{d{x}^{\left(n\right)}}\left(-2{\tanh }^2\left(x\right)\right)\frac{{𝓉}^n}{n!}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=1\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}{\left(-1\right)}^n\left(\frac{d^{\left(n\right)}}{d{𝓉}^{\left(n\right)}}{\left(-2{\tanh }^2\left(x+𝓉\right)\right)}_{\left|𝓉=0\right.}\right)\frac{{𝓉}^n}{n!}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=1+{\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}\left(\frac{d^{\left(n\right)}}{d{𝓉}^{\left(n\right)}}{\left(-2{\tanh }^2\left(x-𝓉\right)\right)}_{\left|𝓉=0\right.}\right)\frac{x^{𝓉}}{n!}=1-2{\tanh }^2\left(x-𝓉\right),\hfill \end{array}$$

(19)

and which is totally in agreement with the exact solution.

Example 2.

Consider the Burgers’ system of nonlinear time fractional IVP:

$$\begin{array}{c}\hfill \frac{{\partial }^{𝒶}\mathcal{U}}{\partial {𝓉}^{𝒶}}-\frac{{\partial }^2\mathcal{U}}{\partial x^2}-2\mathcal{U}\frac{\partial \mathcal{U}}{\partial x}+\mathcal{U}\frac{\partial \mathcal{V}}{\partial x}+\mathcal{V}\frac{\partial \mathcal{U}}{\partial x}=0,\hfill \\ \hfill \frac{{\partial }^{𝒶}\mathcal{V}}{\partial {𝓉}^{𝒶}}-\frac{{\partial }^2\mathcal{V}}{\partial x^2}-2\mathcal{V}\frac{\partial \mathcal{V}}{\partial x}++\mathcal{U}\frac{\partial \mathcal{V}}{\partial x}+\mathcal{V}\frac{\partial \mathcal{U}}{\partial x}=0,\hfill \\ \mathit{subject} \mathit{to} \mathit{ICs}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\mathcal{U}\left(x,0\right)=\mathit{sin}\left(x\right) \mathrm{and} \mathcal{V}\left(x,0\right)=\mathit{sin}\left(x\right), \hfill \end{array}$$

(20)

where $𝒶\in \left(0,1\right]$and$\left(x,𝓉\right)\in ℝ\times \left[0,1\right]$. The exact solutions when$𝒶=1$, is$\mathcal{U}\left(x,𝓉\right)=\sin \left(x\right)e^{-x}$and$\mathcal{V}\left(x,𝓉\right)=\sin \left(x\right)e^{-x}$.

By taking the Laplace transform operator on both sides of (20) and using the second part of Lemma 2 and the initial conditions of $\left(20\right)$, the Laplace fractional equations will be:

$$\begin{array}{cc}\mathfrak{U}\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}+\hfill & \frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\right\}+\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\bigvee \right\}\right\}\hfill \\ & -\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\bigvee \right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\right\},\hfill \\ \bigvee \left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}+\hfill & \frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\bigvee \right\}\right\}+\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\bigvee \right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\bigvee \right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\bigvee \right\}\right\}\hfill \\ & -\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\bigvee \right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}\right\},\hfill \end{array}$$

(21)

where $\mathfrak{U}\left(x,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{U}\left(x,𝓉\right)\right]$ and $\bigvee \left(x,𝓎,\mathfrak{s}\right)= \mathcal{L}\left[\mathcal{V}\left(x,𝓉\right)\right]$.

According to the last discussion of the proposed method, the $k-th$ Laplace series solutions, ${\mathfrak{U}}_k\left(x,\mathfrak{s}\right)$ and ${\bigvee }_k\left(x,\mathfrak{s}\right)$ for (21) are expressed as:

$$\begin{array}{c}{\mathfrak{U}}_k\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{𝒽}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}},\\ \hspace{1em}{\bigvee }_k\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{ℊ}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}}.\end{array}$$

(22)

As well we define the $k-th$ Laplace residual functions of (21) are:

$$\begin{array}{cc}\hfill \mathcal{L}\left(Re{s}_{{\mathfrak{U}}_k}\left(x,\mathfrak{s}\right)\right)=& {\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{𝒽}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\right\}-\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\right\},\hfill \\ \hspace{1em}\mathcal{L}\left(Re{s}_{{\bigvee }_k}\left(x,\mathfrak{s}\right)\right)=\hfill & {\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{ℊ}_n\left(x\right)}{{\mathfrak{s}}^{n𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\right\}-\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\right\}\\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{x_k\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{{\bigvee }_k\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}\right\}.\hfill \end{array}$$

(23)

By letting $k=1$, in (23), the $1-st$ Laplace residual functions are:

$$\begin{array}{c}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_1}\left(x,\mathfrak{s}\right)\right)=\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}-\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\frac{1}{{\mathfrak{s}}^{𝒶+1}}\left({𝒽}_1\left(x\right)+\sin \left(x\right)\right)+\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left(\cos \left(x\right)\left({ℊ}_1\left(x\right)-{𝒽}_1\left(x\right)\right)+\sin \left(x\right)\left({{ℊ}^{\prime }}_1\left(x\right)-{{𝒽}^{\prime }}_1\left(x\right)\right)-{𝒽}_1^{\prime \prime }\left(x\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left({𝒽}_1\left(x\right){{ℊ}^{\prime }}_1\left(x\right)+{ℊ}_1\left(x\right){{𝒽}^{\prime }}_1\left(x\right)-2{𝒽}_1\left(x\right){{𝒽}^{\prime }}_1\left(x\right)\right)\frac{\mathsf{\Gamma }\left(2𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(𝒶+1\right)\right)}^2},\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \mathcal{L}\left(Re{s}_{x_1}\left(x,\mathfrak{s}\right)\right)=\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} -\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{g_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{ℊ}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}+\frac{{𝒽}_1\left(x\right)}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\frac{1}{{\mathfrak{s}}^{𝒶+1}}\left({ℊ}_1\left(x\right)+\sin \left(x\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left(\cos \left(x\right)\left({𝒽}_1\left(x\right)-{ℊ}_1\left(x\right)\right)+\sin \left(x\right)\left({{𝒽}^{\prime }}_1\left(x\right)-{{ℊ}^{\prime }}_1\left(x\right)\right)-{{ℊ}_1}^{″}\left(x\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left({ℊ}_1\left(x\right){{𝒽}^{\prime }}_1\left(x\right)+{𝒽}_1\left(x\right){{ℊ}^{\prime }}_1\left(x\right)-2{ℊ}_1\left(x\right){{ℊ}^{\prime }}_1\left(x\right)\right)\frac{\mathsf{\Gamma }\left(2𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(𝒶+1\right)\right)}^2}.\hfill \end{array}$$

(24)

To find the $1-st$ Laplace series solution of (21), we simply take the next process $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_1}\left(x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\bigvee }_1}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right)$, which yields that ${𝒽}_1\left(x\right)=-\sin \left(x\right)$ and ${ℊ}_1\left(x\right)=-\sin \left(x\right)$. Hence, the $1-st$ Laplace series solutions of (21) are:

$$\begin{array}{l}{\mathfrak{U}}_1\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}},\\ \hspace{1em}{\bigvee }_1\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}.\end{array}$$

(25)

By letting $k=2$, in (23), the $2-nd$ Laplace residual functions are:

$$\begin{array}{c}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_2}\left(x,\mathfrak{s}\right)\right)=-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} -\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left({𝒽}_2\left(x\right)-\sin \left(x\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left(\cos \left(x\right)\left({ℊ}_2\left(x\right)-{𝒽}_2\left(x\right)\right)+\sin \left(x\right)\left({{ℊ}^{\prime }}_2\left(x\right)-{{𝒽}^{\prime }}_2\left(x\right)\right)-{𝒽}_2^{″}\left(x\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{4𝒶+1}}\left(\cos \left(x\right)\left({𝒽}_2\left(x\right)-{ℊ}_2\left(x\right)\right)+\sin \left(x\right)\left({{𝒽}^{\prime }}_2\left(x\right)-{{ℊ}^{\prime }}_2\left(x\right)\right)\right)\frac{\mathsf{\Gamma }\left(3𝒶+1\right)}{\mathsf{\Gamma }\left(2𝒶+1\right)\mathsf{\Gamma }\left(𝒶+1\right)}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{5𝒶+1}}\left({𝒽}_2\left(x\right){{ℊ}^{\prime }}_2\left(x\right)+{ℊ}_2\left(x\right){{𝒽}^{\prime }}_2\left(x\right)-2{𝒽}_2\left(x\right){{𝒽}^{\prime }}_2\left(x\right)\right)\frac{\mathsf{\Gamma }\left(4𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(2𝒶+1\right)\right)}^2},\hfill \\ \hspace{1em} \mathcal{L}\left(Re{s}_{{\mathrm{V}}_2}\left(x,\mathfrak{s}\right)\right)=-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} -\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2\left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left({ℊ}_2\left(x\right)-\sin \left(x\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left(\cos \left(x\right)\left({𝒽}_2\left(x\right)-{ℊ}_2\left(x\right)\right)+\sin \left(x\right)\left({{𝒽}^{\prime }}_2\left(x\right)-{{ℊ}^{\prime }}_2\left(x\right)\right)-{ℊ}_2^{″}\left(x\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{4𝒶+1}}\left(\cos \left(x\right)\left({ℊ}_2\left(x\right)-{𝒽}_2\left(x\right)\right)+\sin \left(x\right)\left({{ℊ}^{\prime }}_2\left(x\right)-{{𝒽}^{\prime }}_2\left(x\right)\right)\right)\frac{\mathsf{\Gamma }\left(3𝒶+1\right)}{\mathsf{\Gamma }\left(2𝒶+1\right)\mathsf{\Gamma }\left(𝒶+1\right)}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{5𝒶+1}}\left({ℊ}_2\left(x\right){{𝒽}^{\prime }}_2\left(x\right)+{𝒽}_2\left(x\right){{ℊ}^{\prime }}_2\left(x\right)-2{ℊ}_2\left(x\right){{ℊ}^{\prime }}_2\left(x\right)\right)\frac{\mathsf{\Gamma }\left(4𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(2𝒶+1\right)\right)}^2}.\hfill \end{array}$$

(26)

To find the $2-nd$ Laplace series solution of (21), we simply find out the next process $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{2𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_2}\left(x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\mathrm{V}}_2}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right),$ and by solving limits, we obtain ${𝒽}_2\left(x\right)=\sin \left(x\right)$ and ${ℊ}_2\left(x\right)=\sin \left(x\right)$. Hence, the $2-nd$ Laplace series solutions of (21) are:

$$\begin{array}{l}{\mathfrak{U}}_2\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}},\\ \hspace{1em}{\mathrm{\bigvee }}_2\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}.\end{array}$$

(27)

Similarly, for $k=3$, we have:

$$\begin{array}{c}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_3}\left(x,\mathfrak{s}\right)\right)=-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} -\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\},\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\mathcal{L}\left(Re{s}_{{\mathrm{V}}_3}\left(x,\mathfrak{s}\right)\right)=-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{\frac{{\partial }^2}{\partial x^2}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{2}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{ℊ}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\frac{\partial }{\partial x}{\mathcal{L}}^{-1}\left\{\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{𝒽}_3\left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}.\hfill \end{array}$$

(28)

By solving $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{3𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_3}\left(x,\mathfrak{s}\right)\right), \mathcal{L}\left(Re{s}_{{\mathrm{V}}_3}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right)$. It yields that: ${𝒽}_3\left(x\right)=-\sin \left(x\right)$ and ${ℊ}_3\left(x\right)=-\sin \left(x\right)$. Hence, the $3-rd$ Laplace series solutions of (21) are:

$$\begin{array}{l}{\mathfrak{U}}_3\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{3𝒶+1}},\\ \hspace{1em}{\mathrm{\bigvee }}_3\left(x,\mathfrak{s}\right)=\frac{\sin \left(x\right)}{\mathfrak{s}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{\sin \left(x\right)}{{\mathfrak{s}}^{2𝒶+1}}-\frac{\sin \left(x\right)}{{\mathfrak{s}}^{3𝒶+1}}.\end{array}$$

(29)

Using Mathematica, we can process the above steps for any $k$, and by the fact that $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{k𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_k}\left(x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\mathrm{V}}_k}\left(x,\mathfrak{s}\right)\right)\right)=\left(0,0\right)$, one can obtain that ${𝒽}_k\left(x\right)={\left(-1\right)}^k\sin \left(x\right)$ and ${ℊ}_k\left(x\right)={\left(-1\right)}^k\sin \left(x\right).$ Thus, the $k-th$ Laplace series solutions of (21) could be formulated on the fractional expansion:

$$\begin{gathered} {\mathfrak{U}}_k\left(x,\mathfrak{s}\right)=\sin \left(x\right)\left(\frac{1}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶+1}}+\frac{1}{{\mathfrak{s}}^{2𝒶+1}}-\frac{1}{{\mathfrak{s}}^{3𝒶+1}}+\dots +{\left(-1\right)}^k\frac{1}{{\mathfrak{s}}^{k𝒶+1}}\right)=\sin \left(x\right){\displaystyle \sum }_{n=0}^k\frac{{\left(-1\right)}^n}{{\mathfrak{s}}^{n𝒶+1}}, \\ {\mathrm{V}}_k\left(x,\mathfrak{s}\right)=\sin \left(x\right)\left(\frac{1}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶+1}}+\frac{1}{{\mathfrak{s}}^{2𝒶+1}}-\frac{1}{{\mathfrak{s}}^{3𝒶+1}}+\dots +{\left(-1\right)}^k\frac{1}{{\mathfrak{s}}^{k𝒶+1}}\right)=\sin \left(x\right){\displaystyle \sum }_{n=0}^k\frac{{\left(-1\right)}^n}{{\mathfrak{s}}^{n𝒶+1}}. \end{gathered}$$

(30)

In the end, we take the inverse LT for the obtained expansions (30) to get that the $k-th$ approximate solutions of the nonlinear system of time-FPDEs (20) have the form:

$$\begin{array}{l}{\mathcal{U}}_k\left(x,𝓉\right)=\sin \left(x\right){\displaystyle {\displaystyle \sum }_{n=0}^k}\frac{{\left(-1\right)}^n{𝓉}^{n𝒶}}{\mathsf{\Gamma }\left(n𝒶+1\right)},\\ \hspace{1em}{\mathcal{\bigvee }}_k\left(x,𝓉\right)=\sin \left(x\right){\displaystyle {\displaystyle \sum }_{n=0}^k}\frac{{\left(-1\right)}^n{𝓉}^{nx}}{\mathsf{\Gamma }\left(n𝒶+1\right)}.\end{array}$$

(31)

When $k\to \infty$ and $𝒶=1$ in (31), the Maclaurin series expansions of the closed forms are:

$$\begin{array}{l}\mathcal{U}\left(x,t\right)=\sin \left(x\right)e^{-𝓉},\\ \hspace{1em}\mathcal{V}\left(x,𝓉\right)=\sin \left(x\right)e^{-𝓉}.\end{array}$$

(32)

and which is totally in agreement with the exact solution.

Example 3.Consider the Burgers’ system of nonlinear time-FPDEs:

$$\begin{array}{c}\frac{{\partial }^{𝒶}\mathcal{U}}{\partial x^{𝓉}}+\frac{\partial \mathcal{V}}{\partial x}\frac{\partial \mathcal{W}}{\partial 𝓎}-\frac{\partial \mathcal{V}}{\partial 𝓎}\frac{\partial \mathcal{W}}{\partial x}+\mathcal{U}=0,\\ \frac{{\partial }^{𝒶}\mathcal{V}}{\partial {𝓉}^{\beta }}+\frac{\partial \mathcal{U}}{\partial x}\frac{\partial \mathcal{W}}{\partial 𝓎}+\frac{\partial \mathcal{U}}{\partial 𝓎}\frac{\partial \mathcal{W}}{\partial x}-\mathcal{V}=0,\\ \frac{{\partial }^{𝒶}\mathcal{W}}{\partial {𝓉}^{\gamma }}+\frac{\partial \mathcal{U}}{\partial x}\frac{\partial \mathcal{V}}{\partial 𝓎}+\frac{\partial \mathcal{U}}{\partial 𝓎}\frac{\partial \mathcal{V}}{\partial x}-\mathcal{W}=0,\\ \mathit{subject} \mathit{to} \mathit{ICs}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ \hspace{1em}\hspace{1em}\hspace{1em}\mathcal{U}\left(x,𝓎,0\right)=e^{x+𝓎}, \mathcal{V}\left(x,𝓎,0\right)=e^{x-𝓎} \mathrm{and} \mathcal{W}\left(x,𝓎,0\right)=e^{-x+𝓎},\end{array}$$

(33)

where$𝒶\in \left(0,1\right]$and$\left(x,𝓎,𝓉\right)\in {ℝ}^2\times \left[0,1\right]$. The exact solutions when$𝒶=1$, are$\left(\mathcal{U}\left(x,𝓎,𝓉\right),\mathcal{V}\left(x,𝓎,𝓉\right),\mathcal{W}\left(x,𝓎,𝓉\right)\right)=\left(e^{x+𝓎-𝓉},e^{x-𝓎+𝓉},e^{-x+𝓎+𝓉}\right).$

By taking the LT operator on both sides of (33) and using the second part of Lemma 2 and the ICs of $\left(33\right)$, the Laplace fractional equations will be:

$$\begin{gathered} \mathfrak{U}\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{x+𝓎}}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\mathrm{V}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\sout{𝒲}\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\mathrm{V}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\sout{𝒲}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathfrak{U}, \\ \mathrm{V}\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{x-𝓎}}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\sout{𝒲}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\sout{𝒲}\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathrm{V}, \\ \sout{𝒲}\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{-𝓎+y}}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\mathrm{V}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\mathfrak{U}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\mathrm{V}\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\sout{𝒲}, \end{gathered}$$

(34)

where $\mathfrak{U}\left(x,𝓎,\mathfrak{s}\right)=\mathcal{L}\left[\mathcal{U}\left(x,𝓎,𝓉\right)\right]$, $\mathrm{V}\left(x,𝓎,\mathfrak{s}\right)=\mathcal{L}\left[\mathcal{V}\left(x,𝓎,𝓉\right)\right]$ and $\sout{𝒲}\left(x,𝓎,\mathfrak{s}\right)=\mathcal{L}\left[\mathcal{W}\left(x,𝓎,𝓉\right)\right]$.

According to the last discussion of the proposed method, the $k-th$ Laplace series solutions ${\mathfrak{U}}_k\left(x,y,\mathfrak{s}\right)$, ${\mathrm{V}}_k\left(x,𝓎,\mathfrak{s}\right)$ and ${\sout{𝒲}}_k\left(x,𝓎,\mathfrak{s}\right)$ for (34) are expressed as:

$$\begin{array}{l}{\mathfrak{U}}_k\left(x,\mathcal{y},\mathfrak{s}\right)=\frac{e^{x+𝓎}}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{𝒽}_n\left(x,𝓎\right)}{{\mathfrak{s}}^{n𝒶+1}},\\ {\mathrm{\bigvee }}_k\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{x-𝓎}}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{ℊ}_n\left(x,𝓎\right)}{{\mathfrak{s}}^{n𝒶+1}},\\ \hspace{1em}{\sout{𝒲}}_k\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{-𝓎+y}}{\mathfrak{s}}+{\displaystyle {\displaystyle \sum }_{n=1}^k}\frac{{\mathfrak{f}}_n\left(x,𝓎\right)}{{\mathfrak{s}}^{n𝒶+1}}.\end{array}$$

(35)

As well, the $k-th$ Laplace fractional residual functions of (34) are defined as:

$$\begin{gathered} \mathcal{L}\left(Re{s}_{{\mathfrak{U}}_k}\left(x,𝓎,\mathfrak{s}\right)\right)={\displaystyle \sum }_{n=1}^k\frac{{𝒽}_n\left(x,𝓎\right)}{{\mathfrak{s}}^{n𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{{\mathrm{V}}_k\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{{\sout{𝒲}}_k\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{{\mathrm{V}}_k\right\}{D}_x{\mathcal{L}}^{-1}\left\{{\sout{𝒲}}_k\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}{\mathfrak{U}}_k, \\ \mathcal{L}\left(Re{s}_{{\mathrm{V}}_k}\left(x,𝓎,\mathfrak{s}\right)\right)={\displaystyle \sum }_{n=1}^k\frac{{ℊ}_n\left(x,𝓎\right)}{{\mathfrak{s}}^{n𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{{\sout{𝒲}}_k\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}{D}_x{\mathcal{L}}^{-1}\left\{{\sout{𝒲}}_k\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}{\mathrm{V}}_k, \\ \mathcal{L}\left(Re{s}_{{\sout{𝒲}}_k}\left(x,𝓎,\mathfrak{s}\right)\right)={\displaystyle \sum }_{n=1}^k\frac{{\mathfrak{f}}_n\left(x,𝓎\right)}{{\mathfrak{s}}^{n𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{{\mathrm{V}}_k\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{{\mathfrak{U}}_k\right\}{D}_x{\mathcal{L}}^{-1}\left\{{\mathrm{V}}_k\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}{\sout{𝒲}}_k. \end{gathered}$$

(36)

For $k=1$, in (36), the $1-st$ Laplace residual functions are expressed as:

$$\begin{gathered} \begin{array}{cc}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_1}\left(x,𝓎,\mathfrak{s}\right)\right)\hfill & \\ & =\frac{{𝒽}_1}{{\mathfrak{s}}^{𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{{ℊ}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{{\mathfrak{f}}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & -\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{{ℊ}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{{\mathfrak{f}}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}+\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{x+𝓎}}{\mathfrak{s}}+\frac{h_1}{{\mathfrak{s}}^{𝒶+1}}\right)\hfill \\ & =\frac{1}{{\mathfrak{s}}^{𝒶+1}}\left({𝒽}_1+e^{x+𝓎}\right)+\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left[e^{x-𝓎}\left(\frac{\partial f_1}{\partial x}+\frac{\partial {\mathfrak{f}}_1}{\partial 𝓎}\right)+e^{-x+𝓎}\left(\frac{\partial {ℊ}_1}{\partial x}+\frac{\partial {ℊ}_1}{\partial 𝓎}\right)+{𝒽}_1\right]\hfill \\ & +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left(\frac{\partial {ℊ}_1}{\partial x}\frac{\partial {\mathfrak{f}}_1}{\partial 𝓎}-\frac{\partial {\mathfrak{f}}_1}{\partial x}\frac{\partial {ℊ}_1}{\partial 𝓎}\right)\frac{\mathsf{\Gamma }\left(2𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(𝒶+1\right)\right)}^2},\hfill \end{array} \\ \begin{array}{cc}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathcal{L}\left(Re{s}_{{\mathrm{V}}_1}\left(x,𝓎,\mathfrak{s}\right)\right)\hfill & \\ & =\frac{{ℊ}_1}{{\mathfrak{s}}^{𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}+\frac{{𝒽}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{{\mathfrak{f}}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}+\frac{{𝒽}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{{\mathfrak{f}}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{{ℊ}_1}{{\mathfrak{s}}^{𝒶+1}}\right)\hfill \\ & =\frac{1}{{\mathfrak{s}}^{𝒶+1}}\left({ℊ}_1-e^{x-𝓎}\right)+\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left[e^{x+𝓎}\left(\frac{\partial {\mathfrak{f}}_1}{\partial x}+\frac{\partial {\mathfrak{f}}_1}{\partial 𝓎}\right)+e^{-x+𝓎}\left(\frac{\partial {𝒽}_1}{\partial x}-\frac{\partial {𝒽}_1}{\partial 𝓎}\right)-{ℊ}_1\right]\hfill \\ & +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left(\frac{\partial {ℊ}_1}{\partial x}\frac{\partial {\mathfrak{f}}_1}{\partial 𝓎}+\frac{\partial {\mathfrak{f}}_1}{\partial x}\frac{\partial {𝒽}_1}{\partial 𝓎}\right)\frac{\mathsf{\Gamma }\left(2\beta +1\right)}{{\left(\mathsf{\Gamma }\left(\beta +1\right)\right)}^2},\hfill \end{array} \\ \begin{array}{cc}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathcal{L}\left(Re{s}_{{\sout{𝒲}}_1}\left(x,𝓎,\mathfrak{s}\right)\right)\hfill & \\ & =\frac{{\mathfrak{f}}_1\left(x,𝓎\right)}{{\mathfrak{s}}^{𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}+\frac{{𝒽}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{{ℊ}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}+\frac{{𝒽}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{𝓎-y}}{\mathfrak{s}}+\frac{{ℊ}_1}{{\mathfrak{s}}^{𝒶+1}}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{{\mathfrak{f}}_1}{{\mathfrak{s}}^{𝒶+1}}\right)\hfill \\ & =\frac{1}{{\mathfrak{s}}^{𝒶+1}}\left(f_1-e^{-x+𝓎}\right)+\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left[e^{𝓎+y}\left(\frac{\partial {ℊ}_1}{\partial x}+\frac{\partial {ℊ}_1}{\partial 𝓎}\right)-e^{x-𝓎}\left(\frac{\partial {𝒽}_1}{\partial x}-\frac{\partial {𝒽}_1}{\partial 𝓎}\right)+{\mathfrak{f}}_1\right]\hfill \\ & +\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left(\frac{\partial {𝒽}_1}{\partial x}\frac{\partial {ℊ}_1}{\partial 𝓎}+\frac{\partial {ℊ}_1}{\partial x}\frac{\partial {𝒽}_1}{\partial 𝓎}\right)\frac{\mathsf{\Gamma }\left(2\gamma +1\right)}{{\left(\mathsf{\Gamma }\left(\gamma +1\right)\right)}^2}.\hfill \end{array} \end{gathered}$$

(37)

To find the $1-st$ Laplace series solution of (34), we simply take the next process $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_1}\left(x,𝓎,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\mathrm{V}}_1}\left(x,𝓎,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\mathcal{W}}_1}\left(x,𝓎,\mathfrak{s}\right)\right)\right)=\left(0,0,0\right)$, which yields that ${𝒽}_1\left(x,𝓎\right)=-e^{x+𝓎}$, ${ℊ}_1\left(x,𝓎\right)=e^{x-𝓎}$ and ${\mathfrak{f}}_1\left(x,𝓎\right)=e^{-x+𝓎}$. Hence, the $1-st$ Laplace series solutions of (34) are:

$$\begin{gathered} {\mathfrak{U}}_1\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{x+𝓎}}{\mathfrak{s}}-\frac{e^{x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}, \\ {\mathrm{\bigvee }}_1\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{e^{x-𝓎}}{{\mathfrak{s}}^{𝒶+1}}, \\ \hspace{1em}\hspace{1em}{\sout{𝒲}}_1\left(x,𝓎,\mathfrak{s}\right)=\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{e^{-x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}. \end{gathered}$$

(38)

For $k=2$, in (36), the $2-nd$ Laplace residual functions are:

$$\begin{gathered} \begin{array}{cc}\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_2}\left(x,𝓎,\mathfrak{s}\right)\right)\hfill & \\ & =-\frac{e^{x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2}{{\mathfrak{s}}^{2𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{e^{x-𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}{D}_{𝓎\prime }{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{e^{-x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{\mathfrak{f}}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ & -\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{e^{x-𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{e^{-x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{\mathfrak{f}}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{x+𝓎}}{\mathfrak{s}}-\frac{e^{x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2}{{\mathfrak{s}}^{2𝒶+1}}\right)\hfill \\ & =\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left({𝒽}_2-e^{x+𝓎}\right)+\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left[e^{x-y}\left(\frac{\partial {\mathfrak{f}}_2}{\partial x}+\frac{\partial {\mathfrak{f}}_2}{\partial 𝓎}\right)+e^{-x+𝓎}\left(\frac{\partial {ℊ}_2}{\partial x}+\frac{\partial {ℊ}_2}{\partial 𝓎}\right)+{𝒽}_2\right]\hfill \\ & +\frac{1}{{\mathfrak{s}}^{4𝒶+1}}\left[e^{x-𝓎}\left(\frac{\partial {\mathfrak{f}}_2}{\partial x}+\frac{\partial {\mathfrak{f}}_2}{\partial 𝓎}\right)+e^{-x+𝓎}\left(\frac{\partial {ℊ}_2}{\partial x}+\frac{\partial {ℊ}_2}{\partial 𝓎}\right)\right]\frac{\mathsf{\Gamma }\left(3𝒶+1\right)}{\mathsf{\Gamma }\left(2𝒶+1\right)\mathsf{\Gamma }\left(𝒶+1\right)}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{5𝒶+1}}\left(\frac{\partial g_2}{\partial x}\frac{\partial {\mathfrak{f}}_2}{\partial 𝓎}-\frac{\partial {\mathfrak{f}}_2}{\partial x}\frac{\partial {ℊ}_2}{\partial 𝓎}\right)\frac{\mathsf{\Gamma }\left(4𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(2𝒶+1\right)\right)}^2},\hfill \end{array} \\ \begin{array}{cc}\mathcal{L}\left(Re{s}_{{\mathrm{V}}_2}\left(x,𝓎,\mathfrak{s}\right)\right)\hfill & \\ & =\frac{e^{x-𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2}{{\mathfrak{s}}^{2𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}-\frac{e^{x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{e^{-x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{\mathfrak{f}}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}-\frac{e^{x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{e^{-x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{\mathfrak{f}}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{e^{x-𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2}{{\mathfrak{s}}^{2𝒶+1}}\right)\hfill \\ & =\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left({ℊ}_2-e^{x-𝓎}\right)+\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left[e^{x+𝓎}\left(\frac{\partial {\mathfrak{f}}_2}{\partial x}+\frac{\partial {\mathfrak{f}}_2}{\partial 𝓎}\right)+e^{-x+𝓎}\left(\frac{\partial {𝒽}_2}{\partial x}-\frac{\partial {𝒽}_2}{\partial 𝓎}\right)+g_2\right]\hfill \\ & +\frac{1}{{\mathfrak{s}}^{4𝒶+1}}\left[-e^{x+𝓎}\left(\frac{\partial {\mathfrak{f}}_2}{\partial x}+\frac{\partial {\mathfrak{f}}_2}{\partial 𝓎}\right)+e^{-x+𝓎}\left(\frac{\partial {𝒽}_2}{\partial x}-\frac{\partial {𝒽}_2}{\partial 𝓎}\right)\right]\frac{\mathsf{\Gamma }\left(3𝒶+1\right)}{\mathsf{\Gamma }\left(2𝒶+1\right)\mathsf{\Gamma }\left(𝒶+1\right)}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{5𝒶+1}}\left(\frac{\partial 𝒽}{\partial x}\frac{\partial {\mathfrak{f}}_2}{\partial 𝓎}+\frac{\partial {\mathfrak{f}}_2}{\partial x}\frac{\partial {𝒽}_2}{\partial 𝓎}\right)\frac{\mathsf{\Gamma }\left(4𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(2𝒶+1\right)\right)}^2},\hfill \end{array} \\ \begin{array}{cc}\mathcal{L}\left(Re{s}_{{\sout{𝒲}}_2}\left(x,𝓎,\mathfrak{s}\right)\right)\hfill & \\ & =\frac{e^{-x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{\mathfrak{f}}_2}{{\mathfrak{s}}^{2𝒶+1}}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}-\frac{e^{x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}{D}_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{e^{x-𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{𝓎}{\mathcal{L}}^{-1}\left\{\frac{e^{x+𝓎}}{\mathfrak{s}}-\frac{e^{x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{𝒽}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{x-𝓎}}{\mathfrak{s}}+\frac{e^{x-𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{ℊ}_2}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{-x+𝓎}}{\mathfrak{s}}+\frac{e^{-x+𝓎}}{{\mathfrak{s}}^{𝒶+1}}+\frac{{\mathfrak{f}}_2}{{\mathfrak{s}}^{2𝒶+1}}\right)\hfill \\ & =\hfill \\ & =\frac{1}{{\mathfrak{s}}^{2𝒶+1}}\left({\mathfrak{f}}_2-e^{-x+𝓎}\right)+\frac{1}{{\mathfrak{s}}^{3𝒶+1}}\left[e^{x+𝓎}\left(\frac{\partial {ℊ}_2}{\partial x}+\frac{\partial {ℊ}_2}{\partial 𝓎}\right)-e^{x-𝓎}\left(\frac{\partial {𝒽}_2}{\partial x}-\frac{\partial {𝒽}_2}{\partial 𝓎}\right)+{\mathfrak{f}}_2\right]\hfill \\ & -\frac{1}{{\mathfrak{s}}^{4𝒶+1}}\left[e^{x+𝓎}\left(\frac{\partial {ℊ}_2}{\partial x}+\frac{\partial {ℊ}_2}{\partial 𝓎}\right)+e^{x-𝓎}\left(\frac{\partial {𝒽}_2}{\partial x}-\frac{\partial {𝒽}_2}{\partial 𝓎}\right)\right]\frac{\mathsf{\Gamma }\left(3\gamma +1\right)}{\mathsf{\Gamma }\left(2𝒶+1\right)\mathsf{\Gamma }\left(𝒶+1\right)}\hfill \\ & +\frac{1}{{\mathfrak{s}}^{5𝒶+1}}\left(\frac{\partial {ℊ}_2}{\partial x}\frac{\partial {𝒽}_2}{\partial 𝓎}+\frac{\partial {𝒽}_2}{\partial x}\frac{\partial {ℊ}_2}{\partial 𝓎}\right)\frac{\mathsf{\Gamma }\left(4𝒶+1\right)}{{\left(\mathsf{\Gamma }\left(2𝒶+1\right)\right)}^2}.\hfill \end{array} \end{gathered}$$

(39)

To find the $2-nd$ Laplace series solution of (34), we simply find out the next process $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{2𝒶+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_2}\left(x,𝓎,\mathfrak{s}\right)\right), \mathcal{L}\left(Re{s}_{{\mathrm{V}}_2}\left(x,𝓎,\mathfrak{s}\right)\right), \mathcal{L}\left(Re{s}_{{\mathcal{W}}_2}\left(x,𝓎,\mathfrak{s}\right)\right)\right)=\left(0,0,0\right),$ and by solving limits, we get ${𝒽}_2\left(x,𝓎\right)=e^{x+𝓎}$, ${ℊ}_2\left(x,𝓎\right)=e^{x-𝓎}$ and ${\mathfrak{f}}_2\left(x,𝓎\right)=e^{-x+𝓎}$. Hence, the $2-nd$ Laplace series solutions of (34) are:

$$\begin{array}{c}{\mathfrak{U}}_2\left(x,\mathcal{y},\mathfrak{s}\right)=\frac{e^{x+\mathcal{y}}}{\mathfrak{s}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}},\\ {\bigvee }_2\left(x,\mathcal{y},\mathfrak{s}\right)=\frac{e^{x-\mathcal{y}}}{\mathfrak{s}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}},\\ \hspace{1em}\hspace{1em}\hspace{1em}{\sout{𝒲}}_2\left(x,\mathcal{y},\mathfrak{s}\right)=\frac{e^{-x+\mathcal{y}}}{\mathfrak{s}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}.\end{array}$$

(40)

Similarly, for $k=3$, we have:

$$\begin{array}{l}\mathcal{L}\left(R,e,s_{{\mathfrak{U}}_3},\left(x,\mathcal{y},\mathfrak{s}\right)\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{-e^{x+x}}{{\mathfrak{s}}^{x+1}}+\frac{e^{x+x}}{{\mathfrak{s}}^{2x+1}}+\frac{x_3}{{\mathfrak{s}}^{3x+1}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x-\mathcal{y}}}{\mathfrak{s}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{g}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}{D}_{\mathcal{y}}{\mathcal{L}}^{-1}\left\{\frac{e^{-x+\mathcal{y}}}{\mathfrak{s}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathfrak{f}}_3}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{\mathcal{y}}{\mathcal{L}}^{-1}\left\{\frac{e^{x-\mathcal{y}}}{\mathfrak{s}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{g}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{-x+\mathcal{y}}}{\mathfrak{s}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathfrak{f}}_3}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{x+\mathcal{y}}}{\mathfrak{s}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{h}}_3}{{\mathfrak{s}}^{\mathrm{a}x+1}}\right),\\ \hspace{1em}\hspace{1em}\mathcal{L}\left(Re{s}_{{\bigvee }_3}\left(x,\mathcal{y},\mathfrak{s}\right)\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{g}}_3}{{\mathfrak{s}}^{3𝒶+1}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{\mathcal{y}}{\mathcal{L}}^{-1}\left\{\frac{e^{x+\mathcal{y}}}{\mathfrak{s}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{h}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}{D}_{\mathcal{y}}{\mathcal{L}}^{-1}\left\{\frac{e^{-x+\mathcal{y}}}{\mathfrak{s}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathfrak{f}}_3}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{\mathcal{y}}{\mathcal{L}}^{-1}\left\{\frac{e^{x+\mathcal{y}}}{\mathfrak{s}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{h}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{-x+\mathcal{y}}}{\mathfrak{s}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathfrak{f}}_3}{{\mathfrak{s}}^{2𝒶+1}}\right\}\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{x-\mathcal{y}}}{\mathfrak{s}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{g}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right),\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathcal{L}\left(Re{s}_{{\sout{𝒲}}_3}\left(x,\mathcal{y},\mathfrak{s}\right)\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathfrak{f}}_3}{{\mathfrak{s}}^{2𝒶+1}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_x{\mathcal{L}}^{-1}\left\{\frac{e^{x+\mathcal{y}}}{\mathfrak{s}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{h}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}{D}_{\mathcal{y}}{\mathcal{L}}^{-1}\left\{\frac{e^{x-\mathcal{y}}}{\mathfrak{s}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{g}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{{\mathfrak{s}}^{𝒶}}\mathcal{L}\left\{D_{\mathcal{y}}{\mathcal{L}}^{-1}\left\{\frac{e^{x+\mathcal{y}}}{\mathfrak{s}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{h}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}{D}_x{\mathcal{L}}^{-1}\left\{\frac{e^{x-\mathcal{y}}}{\mathfrak{s}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathcal{g}}_3}{{\mathfrak{s}}^{3𝒶+1}}\right\}\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{1}{{\mathfrak{s}}^{𝒶}}\left(\frac{e^{-x+\mathcal{y}}}{\mathfrak{s}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{{\mathfrak{f}}_3}{{\mathfrak{s}}^{2𝒶+1}}\right).\end{array}$$

(41)

By solving $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{3x+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_3}\left(x,x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{x_3}\left(x,x,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\mathcal{W}}_3}\left(x,x,\mathfrak{s}\right)\right)\right)=\left(0,0,0\right)$, it yields that: $x_3\left(x,y\right)=-e^{x+y}$, $x_3\left(x,x\right)=e^{x-x}$ and ${\mathfrak{f}}_3\left(x,x\right)=e^{-x+x}$. Hence, the $3-rd$ Laplace series solutions of (34) are:

$$\begin{array}{c}{\mathfrak{U}}_3\left(x,\mathcal{y},\mathfrak{s}\right)=\frac{e^{x+\mathcal{y}}}{\mathfrak{s}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}-\frac{e^{x+\mathcal{y}}}{{\mathfrak{s}}^{3𝒶+1}},\\ {\bigvee }_3\left(x,\mathcal{y},\mathfrak{s}\right)=\frac{e^{x-\mathcal{y}}}{\mathfrak{s}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{e^{x-\mathcal{y}}}{{\mathfrak{s}}^{3𝒶+1}},\\ \hspace{1em}\hspace{1em}\hspace{1em}{\sout{𝒲}}_3\left(x,\mathcal{y},\mathfrak{s}\right)=\frac{e^{-x+\mathcal{y}}}{\mathfrak{s}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{2𝒶+1}}+\frac{e^{-x+\mathcal{y}}}{{\mathfrak{s}}^{3𝒶+1}}.\end{array}$$

(42)

Using Mathematica, we can process the above steps for any $k$, and by the fact that $\underset{\mathfrak{s}\to \infty }{\lim }{\mathfrak{s}}^{k𝓪+1}\left(\mathcal{L}\left(Re{s}_{{\mathfrak{U}}_k}\left(𝓍,𝓎,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\vee }_k}\left(𝓍,𝓎,\mathfrak{s}\right)\right),\mathcal{L}\left(Re{s}_{{\mathcal{W}}_k}\left(𝓍,𝓎,\mathfrak{s}\right)\right)\right)=\left(0,0,0\right)$, one can obtain that ${𝒽}_k\left(𝓍,𝓎\right)={\left(-1\right)}^k{e}^{𝓍+𝓎},$ ${𝓰}_k\left(𝓍,𝓎\right)=e^{𝓍-𝓎}$ and ${\mathfrak{f}}_k\left(𝓍,𝓎\right)=e^{-𝓍+𝓎}.$ Thus, the $k$th-Laplace series solutions of (34) could be formulated by the fractional expansions:

$$\begin{gathered} {\mathfrak{U}}_k\left(𝓍,𝓎,\mathfrak{s}\right)=e^{𝓍+𝓎}\left(\frac{1}{\mathfrak{s}}-\frac{1}{{\mathfrak{s}}^{𝓪+1}}+\frac{1}{{\mathfrak{s}}^{2𝓪+1}}-\frac{1}{{\mathfrak{s}}^{3𝓪+1}}+\dots +{\left(-1\right)}^k\frac{1}{{\mathfrak{s}}^{k𝓪+1}}\right)=e^{𝓍+𝓎}{\displaystyle \sum }_{n=0}^k\frac{{\left(-1\right)}^n}{{\mathfrak{s}}^{n𝓪+1}}, \\ {\vee }_k\left(𝓍,𝓎,\mathfrak{s}\right)=e^{𝓍-𝓎}\left(\frac{1}{\mathfrak{s}}+\frac{1}{{\mathfrak{s}}^{𝓪+1}}+\frac{1}{{\mathfrak{s}}^{2𝓪+1}}+\frac{1}{{\mathfrak{s}}^{3𝓪+1}}+\dots +\frac{1}{{\mathfrak{s}}^{k𝓪+1}}\right)=e^{𝓍-𝓎}{\displaystyle \sum }_{n=0}^k\frac{1}{{\mathfrak{s}}^{n𝓪+1}}, \\ \hspace{1em}\hspace{1em}\hspace{1em}{\sout{𝒲}}_k\left(𝓍,𝓎,\mathfrak{s}\right)=e^{-𝓍+𝓎}\left(\frac{1}{\mathfrak{s}}+\frac{1}{{\mathfrak{s}}^{𝓪+1}}+\frac{1}{{\mathfrak{s}}^{2𝓪+1}}+\frac{1}{{\mathfrak{s}}^{3𝓪+1}}+\dots +\frac{1}{{\mathfrak{s}}^{k𝓪+1}}\right)=e^{-𝓍+𝓎}{\displaystyle \sum }_{n=0}^k\frac{1}{{\mathfrak{s}}^{n𝓪+1}}. \end{gathered}$$

(43)

In the end, we take the inverse LT for the obtained expansions (43) to conclude that the $k-th$ approximate solutions of the nonlinear systems of time-FPDEs (33) have the form:

$$\begin{gathered} {\mathcal{U}}_k\left(𝓍,𝓎,𝓉\right)=e^{𝓍+𝓎}{\displaystyle \sum }_{n=0}^k\frac{{\left(-1\right)}^n{𝓉}^{n𝓪}}{\mathsf{\Gamma }\left(n𝓪+1\right)}, \\ {\mathcal{V}}_k\left(𝓍,𝓎,𝓉\right)=e^{𝓍-𝓎}{\displaystyle \sum }_{n=0}^k\frac{{𝓉}^{n𝓪}}{\mathsf{\Gamma }\left(n𝓪+1\right)}, \\ \hspace{1em}\hspace{1em}\hspace{1em}{\mathcal{W}}_k\left(𝓍,𝓎,𝓉\right)=e^{-𝓍+𝓎}{\displaystyle \sum }_{n=0}^k\frac{{𝓉}^{n𝓪}}{\mathsf{\Gamma }\left(n𝓪+1\right)}. \end{gathered}$$

(44)

When $k\to \infty$ and $𝓪=1$ in (44), the Maclaurin series expansions of the closed forms are:

$$\begin{gathered} \mathcal{U}\left(𝓍,𝓎,𝓉\right)=e^{𝓍+𝓎-𝓉}, \\ \mathcal{V}\left(𝓍,𝓎,𝓉\right)=e^{𝓍-𝓎+𝓉}, \\ \hspace{1em}\hspace{1em}\hspace{1em}\mathcal{W}\left(𝓍,𝓎,𝓉\right)=e^{-𝓍+𝓎+𝓉}, \end{gathered}$$

(45)

and which is totally in agreement with the exact solution.

5. Graphical and Numerical Results

This section deals with the validity and efficiency of the Laplace RPSM for systems of time-FPDEs discussed in Examples 1–3 through different graphical representations and tabulated data for the obtained approximation and exact solutions.

The absolute error functions calculated demonstrate the accuracy of the Laplace RPSM.

Table 1. Numerical results for Example 1 at $𝓍=1, 𝓪=1, \mathrm{and} n=7$.

$\mathbf{Numerical} \mathbf{Results} \mathbf{of} \mathcal{U}\left(\mathit{𝓍},\mathit{𝓉}\right)$
${\mathit{𝓉}}_{\mathit{i}}$	$\mathcal{U}\left(\mathit{𝓍},\mathit{𝓉}\right)$	${\mathcal{U}}_7\left(\mathit{𝓍},\mathit{𝓉}\right)$	$\left|\mathcal{U}-{\mathcal{U}}_7\right|$
$0.1$	$2.4325957403980487$	$2.4325957405419616$	$1.439128816116408\times {10}^{-10}$
$0.2$	$2.3280735405356980$	$2.3280735789915883$	$3.845589047202225\times {10}^{-8}$
$0.3$	$2.2087355542343268$	$2.2087365789980504$	$1.024763723656008\times {10}^{-6}$
$0.4$	$2.0740991339960706$	$2.0741097277878920$	$1.059379182155595\times {10}^{-5}$
$0.5$	$1.9242343145200196$	$1.9242993140796738$	$6.499955965422188\times {10}^{-5}$
$0.6$	$1.7598979245104498$	$1.7601838692172310$	$2.859447067811161\times {10}^{-4}$
$0.7$	$1.5826252249031816$	$1.5836225257595018$	$9.973008563202157\times {10}^{-4}$
$0.8$	$1.3947506404498080$	$1.3976783358533416$	$2.927695403533548\times {10}^{-3}$
$0.9$	$1.1993359892499116$	$1.2068540917157446$	$7.518102465833065\times {10}^{-3}$
Numerical Results of$\mathcal{V}\left(\mathit{𝓍},\mathit{𝓉}\right)$
${\mathit{𝓉}}_{\mathit{i}}$	$\mathcal{V}\left(\mathit{𝓍},\mathit{𝓉}\right)$	${\mathcal{V}}_7\left(\mathit{𝓍},\mathit{𝓉}\right)$	$\left|\mathcal{V}-{\mathcal{V}}_7\right|$
$0.1$	$-0.0261652777033167$	$-0.0261652782859831$	$5.82666359605355\times {10}^{-10}$
$0.2$	$0.11811033546448813$	$0.11811018874091550$	$1.46723572624907\times {10}^{-7}$
$0.3$	$0.26947917996491720$	$0.26947552093723215$	$3.659027685065652\times {10}^{-6}$
$0.4$	$0.42315552517444555$	$0.42312041684151236$	$3.510833293318694\times {10}^{-5}$
$0.5$	$0.57289546593185480$	$0.57269750532341330$	$1.497960608441461\times {10}^{-4}$
$0.6$	$0.71127757216235540$	$0.71048682802189810$	$7.907441404573223\times {10}^{-4}$
$0.7$	$0.83027392365325850$	$0.82780654822770770$	$2.467375425550755\times {10}^{-3}$
$0.8$	$0.92208596593223320$	$0.91572529863721930$	$6.360667295013944\times {10}^{-3}$
$0.9$	$0.98013258169487960$	$0.96613158040482160$	$1.400100129005799\times {10}^{-3}$

,

Table 2. Numerical results for Example 2 at $x=10, 𝓪=1, \mathrm{and} n=7$.

$\mathbf{Numerical} \mathbf{Results} \mathbf{of} \mathcal{U}\left(\mathit{𝓍},\mathit{𝓉}\right)$
${\mathit{𝓉}}_{\mathit{i}}$	$\mathcal{U}\left(\mathit{𝓍},\mathit{𝓉}\right)$	${\mathcal{U}}_7\left(\mathit{𝓍},\mathit{𝓉}\right)$	Abs. Error
$0.2$	$-0.4454068138087739$	$-0.4454068137749855$	$3.378847202029078\times {10}^{-11}$
$0.4$	$-0.3646682560957220$	$-0.3646682476310534$	$8.464668632690575\times {10}^{-9}$
$0.6$	$-0.2985651159368848$	$-0.2985649035620645$	$2.123748202298436\times {10}^{-7}$
$0.8$	$-0.2444444422138206$	$-0.2444423647493085$	$2.077464512112437\times {10}^{-6}$
$1$	$-0.2001341822594486$	$-0.2001220515057324$	$1.213075371614413\times {10}^{-5}$
Numerical Results of$\mathcal{V}\left(\mathit{𝓍},\mathit{𝓉}\right)$
${\mathit{𝓉}}_i$	$\mathcal{V}\left(\mathit{𝓍},\mathit{𝓉}\right)$	${\mathcal{V}}_7\left(\mathit{𝓍},\mathit{𝓉}\right)$	Abs. Error
$0.2$	$-0.4454068138087739$	$-0.4454068137749855$	$3.378847202029078\times {10}^{-11}$
$0.4$	$-0.3646682560957220$	$-0.3646682476310534$	$8.464668632690575\times {10}^{-9}$
$0.6$	$-0.2985651159368848$	$-0.2985649035620645$	$2.123748202298436\times {10}^{-7}$
$0.8$	$-0.2444444422138206$	$-0.2444423647493085$	$2.077464512112437\times {10}^{-6}$
$1$	$-0.2001341822594486$	$-0.2001220515057324$	$1.213075371614413\times {10}^{-5}$

and

Table 3. Numerical results for Example 3 at different values of $𝓍,𝓎=0.4,𝓪=1, \mathrm{and} n=7$.

$\mathbf{Numerical} \mathbf{Results} \mathbf{of} \mathcal{U}\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$
${\mathit{𝓍}}_{\mathit{i}}$	${\mathit{𝓉}}_{\mathit{i}}$	$\mathcal{U}\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$	${\mathcal{U}}_7\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$	$\left|\mathcal{U}-{\mathcal{U}}_7\right|$
$0.5$	$0.3$	$1.8221188003905090$	$1.8221187965176793$	$3.872829834605795\times {10}^{-9}$
	$0.6$	$1.3498588075760032$	$1.3498578473967784$	$9.601792247959650\times {10}^{-7}$
	$0.9$	$1.0$	$0.9999761481976381$	$2.385180236186279\times {10}^{-5}$
$1$	$0.3$	$3.0041660239464334$	$3.0041660175612160$	$6.385217243831676\times {10}^{-9}$
	$0.6$	$2.2255409284924680$	$2.2255393454245560$	$1.583067911870017\times {10}^{-6}$
	$0.9$	$1.6487212707001284$	$1.6486819457262294$	$3.932497389902423\times {10}^{-5}$
$1.5$	$0.3$	$4.9530324243951150$	$4.9530324138676710$	$1.052744380558579\times {10}^{-8}$
	$0.6$	$3.6692966676192444$	$3.6692940575815050$	$2.610037739270154\times {10}^{-6}$
	$0.9$	$2.718281828459045$	$2.718216992538108$	$6.483592093697865\times {10}^{-5}$
Numerical Results of$\mathcal{V}\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$
${\mathit{𝓍}}_i$	${\mathit{𝓉}}_i$	$\mathcal{V}\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$	${\mathcal{V}}_7\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$	$\left|\mathcal{V}-{\mathcal{V}}_7\right|$
$0.5$	$0.3$	$1.4918246976412703$	$1.4918246957811063$	$1.860164022815524\times {10}^{-9}$
	$0.6$	$2.0137527074704766$	$2.0137522144484308$	$4.930220458554402\times {10}^{-7}$
	$0.9$	$2.7182818284590450$	$2.7182687338653180$	$1.309459372711430\times {10}^{-5}$
$1$	$0.3$	$2.4596031111569500$	$2.4596031080900580$	$3.066891629543988\times {10}^{-9}$
	$0.6$	$3.3201169227365480$	$3.3201161098806145$	$8.128559336739727\times {10}^{-7}$
	$0.9$	$4.4816890703380645$	$4.4816674810028550$	$2.158933520934880\times {10}^{-5}$
$1.5$	$0.3$	$4.0551999668446745$	$4.0551999617882260$	$5.056448593165896\times {10}^{-9}$
	$0.6$	$5.4739473917272010$	$5.4739460515543330$	$1.340172867791977\times {10}^{-6}$
	$0.9$	$7.3890560989306500$	$7.3890205041344710$	$3.559479617898376\times {10}^{-5}$
Numerical Results of$\mathcal{W}\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$
${\mathit{𝓍}}_i$	${\mathit{𝓉}}_i$	$\mathcal{W}\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$	${\mathcal{W}}_7\left(\mathit{𝓍},\mathit{𝓎},\mathit{𝓉}\right)$	$\left|\mathcal{W}-{\mathcal{W}}_7\right|$
$0.5$	$0.3$	$1.2214027581601699$	$1.2214027566371966$	$1.52297330302531\times {10}^{-9}$
	$0.6$	$1.6487212707001282$	$1.6487208670478177$	$4.03652310465574\times {10}^{-7}$
	$0.9$	$2.2255409284924674$	$2.2255302075458845$	$1.07209465829427\times {10}^{-5}$
$1$	$0.3$	$0.7408182206817179$	$0.7408182197579880$	$9.2372987037236\times {10}^{-10}$
	$0.6$	$1.0$	$0.9999997551724980$	$2.44827502049461\times {10}^{-7}$
	$0.9$	$1.3498588075760030$	$1.3498523049931992$	$6.50258280376903\times {10}^{-6}$
$1.5$	$0.3$	$0.4493289641172216$	$0.4493289635569510$	$5.60270607685710\times {10}^{-10}$
	$0.6$	$0.6065306597126334$	$0.6065305112172470$	$1.48495386431690\times {10}^{-7}$
	$0.9$	$0.8187307530779818$	$0.8187268090621439$	$3.94401583792003\times {10}^{-6}$

illustrate several values of the approximate and exact solutions as well as the absolute errors for systems of time-FPDEs (7), (20), and (44) at selected grid points in the domain. From the tables, the approximate solutions are harmonic with the exact solutions, which confirms the performance and accuracy of the Laplace RPSM, whilst the accuracy is in advance by using only a few of the Laplace RPS iterations. Further, numerical simulations for the attained results of the problems studied are achieved at various values of $𝓪$ as illustrated in

Table 4. Numerical results of approximated solutions, at $n = 7,$ 𝓍 = 1, and different values of 𝒶, for Example 1.

${\mathcal{U}}_7$(𝓍,𝓉)
𝓉i	𝒶 = 1	𝒶 = 0.97	𝒶 = 0.87	𝒶 = 0.77
0.15	2.3821389434357910	2.3660488221041770	2.3164985673552834	2.2494580999988480
0.30	2.2087365789980504	2.1798840419417136	2.0946691478570347	1.9890784760486064
0.45	2.0010680201298103	1.9605498108679509	1.8475490983553642	1.7270938997910783
0.60	1.7601838692172317	1.7110610843964820	1.5871611481074164	1.4930439371551603
0.75	1.4915784441913985	1.4404546365696143	1.3361912466819796	1.3308640595516880
0.90	1.2068540917157446	1.1653285241400642	1.1284934441148844	1.2968610082469043
${\mathcal{V}}_7$(𝓍,𝓉)
𝓉i	𝒶 = 1	𝒶 = 0.97	𝒶 = 0.87	𝒶 = 0.77
0.15	0.0448459609008298	0.0663869988097889	0.129641467880782	0.2067709053195330
0.30	0.2694755209372319	0.3009312650284730	0.3825827826197465	0.4565703174103681
0.45	0.4988712398275166	0.5298370150916816	0.5914524244945669	0.6089345684206415
0.60	0.7104868280218981	0.7269468972314479	0.7268268374287052	0.6726093727343672
0.75	0.8759824453139846	0.8623611282911388	0.7744858455437431	0.7307371590702421
0.90	0.9661315804048216	0.9112450135222343	0.7596624176539246	0.9841962235906878

,

Table 5. Numerical results of approximated solutions, at n = 7, 𝓍 = 10, and different values of 𝒶, for Example 2.

${\mathcal{U}}_7$(𝓍,𝓉)
𝓉i	𝒶 = 1	𝒶 = 0.95	𝒶 = 0.85	𝒶 = 0.75
0.2	−0.4454068138087739	−0.4366187365008558	−0.41802625689971980	−0.3985414981473607
0.4	−0.3646682560957220	−0.3565811645027616	−0.34129966933185546	−0.3275638100454310
0.6	−0.2985651159368848	−0.2934671682315440	−0.28497017140893940	−0.2786635555415501
0.8	−0.2444444422138206	−0.2429390773523567	−0.24160764152557102	−0.2420836601179408
1.0	−0.2001341822594486	−0.2021160815904085	−0.20726098237141660	−0.2132493814323069
${\mathcal{V}}_7$(𝓍,𝓉)
𝓉i	𝒶 = 1	𝒶 = 0.95	𝒶 = 0.85	𝒶 = 0.75
0.2	−0.4454068138087739	−0.4366187365008558	−0.41802625689971980	−0.3985414981473607
0.4	−0.3646682560957220	−0.3565811645027616	−0.34129966933185546	−0.3275638100454310
0.6	−0.2985651159368848	−0.2934671682315440	−0.28497017140893940	−0.2786635555415501
0.8	−0.2444444422138206	−0.2429390773523567	−0.24160764152557102	−0.2420836601179408
1.0	−0.2001341822594486	−0.2021160815904085	−0.20726098237141660	−0.2132493814323069

and

Table 6. Numerical results of approximated solutions, at y = 0.4 and different values of 𝓍, 𝓉, and 𝒶, with $n = 7$ for Example 3.

${\mathcal{U}}_7$(𝓍,𝓎,𝓉)
𝓍i	𝓉i	𝒶 = 1	𝒶 = 0.95	𝒶 = 0.75	𝒶 = 0.55
0.5	0.3	1.82211879651767	1.781854354922940	1.624165815324987	1.485973706352699
	0.6	1.34985784739677	1.3268102019509882	1.259880792227996	1.224207523196436
	0.9	0.99997614819763	1.0012559086702593	1.025844721313154	1.048930393543767
1	0.3	3.00416601756121	2.937781176251118	2.677796726870322	2.449956457364802
	0.6	2.22553934542455	2.187540202138526	2.077192260692827	2.018376983245085
	0.9	1.64868194572622	1.650791914038841	1.691332012464442	1.729393851319466
1.5	0.3	4.95303241386767	4.8435823139476595	4.414940422202281	4.039295323546481
	0.6	3.66929405758150	3.6066440617774465	3.424711063537949	3.327741064567730
	0.9	2.71821699253810	2.7216957421756143	2.788535064766180	2.851288428088421
${\mathcal{V}}_7$(𝓍,𝓎,𝓉)
𝓍i	𝓉i	𝒶 = 1	𝒶 = 0.95	𝒶 = 0.75	𝒶 = 0.55
0.5	0.3	1.49182469578111	1.533776716329884	1.763906348007679	2.173624234307866
	0.6	2.01375221444843	2.090517467011799	2.500822166419968	3.2070987754409066
	0.9	2.71826873386532	2.836553335181282	3.463909478191494	4.5195588991940581
1	0.3	2.45960310809005	2.528770296717676	2.908189915483245	3.5837005096126577
	0.6	3.32011610988061	3.446680614632507	4.123158700014978	5.2876119683057565
	0.9	4.48166748100285	4.676685819188771	5.711021236474098	7.4514928912833001
1.5	0.3	4.05519996178822	4.169237376913107	4.794794572892834	5.9085232580172781
	0.6	5.47394605155433	5.682615642654406	6.797939451186984	8.7177983233542732
	0.9	7.38902050413447	7.710551386478182	9.415882189994994	12.285434828329576
${\mathcal{W}}_7$(𝓍,𝓎,𝓉)
𝓍i	𝓉i	𝒶 = 1	𝒶 = 0.95	𝒶 = 0.75	𝒶 = 0.55
0.5	0.3	1.22140275663719	1.255750166014239	1.444164372663361	1.7796130062634311
	0.6	1.64872086704781	1.711570940089245	2.047500015627131	2.6257503956122070
	0.9	2.22553020754588	2.322373448258833	2.836009215673681	3.7003018611174463
1	0.3	0.74081821975798	0.761650976626866	0.8759299696849896	1.0793898507221418
	0.6	0.99999975517249	1.038120251437302	1.2418715352399512	1.5925981196913803
	0.9	1.34985230499319	1.408590699671533	1.7201265405336663	2.2443465289594501
1.5	0.3	0.44932896355695	0.46196466932427	0.5312783823751036	0.6546830382456215
	0.6	0.60653051121724	0.62965176096531	0.7532331615474284	0.9659595881935124
	0.9	0.81872680906215	0.85435344633685	1.0433094854190943	1.3612649808335340

.

Numerical comparisons are established to confirm the mathematical results for the obtained approximate solutions supported by the results of numerical comparisons.

Table 7. Numerical comparisons for Example 1 at 𝒶 = 1 and different values of 𝓍 and 𝓉.

$\left|\mathcal{U}-{\mathcal{U}}_2\right|$
	𝓍 = −1		𝓍 = 0.5		𝓍 = 1
𝓉i	LRPSM	MGMLFM [40]	LRPSM	MGMLFM [40]	LRPSM	MGMLFM [40]
0.003	5.59911 × 10−9	4.51907 × 10−5	5.11357 × 10−9	4.53038 × 10−5	5.59013 × 10−9	3.28798 × 10−5
0.006	4.48282 × 10−8	1.79812 × 10−4	4.11223 × 10−8	1.80805 × 10−4	4.46846 × 10−8	1.2751 × 10−4
0.009	1.51413 × 10−7	4.02448 × 10−4	1.3951 × 10−7	4.05897 × 10−4	1.50686 × 10−7	2.78165 × 10−4
$\left|\mathcal{V}-{\mathcal{V}}_2\right|$
	𝓍 = −1		𝓍 = 0.5		𝓍 = 1
𝓉i	LRPSM	MGMLFM [40]	LRPSM	MGMLFM [40]	LRPSM	MGMLFM [40]
0.003	5.94837 × 10−9	1.3969 × 10−5	3.55919 × 10−8	2.14755 × 10−5	6.02338 × 10−9	1.25863 × 10−4
0.006	4.72882 × 10−8	5.5459 × 10−5	2.84908 × 10−7	8.62386 × 10−5	4.84885 × 10−8	4.93352 × 10−4
0.009	1.58592 × 10−7	1.23872 × 10−4	9.62139 × 10−7	1.94798 × 10−4	1.64669 × 10−7	1.08802 × 10−3

shows the absolute errors of the obtained approximate solutions for the system of time-FPDEs (7) at $𝓪=1$ with the absolute errors of the approximate solutions generated by the MGMLFM [40], while

Table 8. Numerical comparisons for Example 1 at different values of 𝓍, 𝒶, and 𝓉 for the function $\mathcal{U}$.

𝒶 = 1
(𝓍,𝓉i)	$\mathcal{U}$(𝓍,𝓉i)	LRPSM	MGMLFM [40]
(−0.5, 0.003)	0.0710535459	0.0710535408	0.0710849
(−0.5, 0.006)	0.0663545199	0.0663544797	0.0664763
(0.5, 0.003)	1.9195090915	1.91950908636	1.91955
(0.5, 0.006)	1.9147708158	1.91477077469	1.91495
𝒶 = 0.9
(𝓍,𝓎,𝓉i)	$\mathcal{U}$(𝓍,𝓎,𝓉i)	LRPSM	MGMLFM [40]
(−0.5, 0.003)	0.0710535459	0.06718985810	0.0671357
(−0.5, 0.006)	0.0663545199	0.05989588986	0.0598679
(0.5, 0.003)	1.9195090915	1.91583860783	1.91561
(0.5, 0.006)	1.9147708158	1.90988083747	1.90838
𝒶 = 0.8
(𝓍,𝓎,𝓉i)	$\mathcal{U}$(𝓍,𝓎,𝓉i)	LRPSM	MGMLFM [40]
(−0.5, 0.003)	0.0710535459	0.06066868120	0.0600753
(−0.5, 0.006)	0.0663545199	0.04976089352	0.0489784
(0.5, 0.003)	1.9195090915	1.90995039533	1.90859
(0.5, 0.006)	1.9147708158	1.89876300851	1.8977

,

Table 9. Numerical comparisons for Example 2 at different values of 𝓍, 𝒶, and 𝓉 for the function $\mathcal{U}$.

𝒶 = 1
(𝓍,𝓉i)	$\mathcal{U}$(𝓍,𝓉i)	LRPSM	MGMLFM [40]
(−10, 0.2)	0.4454068138	0.4454068137	0.4454068
(−10, 0.4)	0.3646682570	0.3646682476	0.3646684
(−5, 0.2)	0.7851007935	0.7851007935	0.7851008
(−5, 0.4)	0.6427861639	0.6427861490	0.6427865
𝒶 = 0.9
(𝓍,𝓎,𝓉i)	$\mathcal{U}$(𝓍,𝓎,𝓉i)	LRPSM	MGMLFM [40]
(−10, 0.2)	0.4454068138	0.42747254328	0.4274714
(−10, 0.4)	0.3646682570	0.34877861288	0.3487726
(−5, 0.2)	0.7851007935	0.75349665629	0.7534867
(−5, 0.4)	0.6427861639	0.61477525733	0.6147676
𝒶 = 0.75
(𝓍,𝓎,𝓉i)	$\mathcal{U}$(𝓍,𝓎,𝓉i)	LRPSM	MGMLFM [40]
(−10, 0.2)	0.4454068138	0.39855149815	0.3985421
(−10, 0.4)	0.3646682570	0.32766381008	0.3275878
(−5, 0.2)	0.7851007935	0.70250317423	0.7024943
(−5, 0.4)	0.6427861639	0.57749363947	0.5774259

and

Table 10. Numerical comparisons for Example 3 at different values of 𝓍, 𝒶, and 𝓉 for the function $\mathcal{U}$.

𝒶 = 1
(𝓍,𝓎,𝓉i)	$\mathcal{U}$(𝓍,𝓎,𝓉i)	LRPSM	MGMLFM [40]	FNDM [41]
(0.5, 0.4, 0.3)	1.8221188	1.8221188	1.8221189	1.82217
(0.5, 0.4, 0.6)	1.3498588	1.3498578	1.3498715	1.35131
(0.5, 0.4, 0.9)	1.0000000	0.9999761	1.00021	1.0105
(1, 0.4, 0.3)	3.0041660	3.0041660	3.0041662	3.00424
(1, 0.4, 0.6)	2.2255400	2.2255394	2.225562	2.22793
(1, 0.4, 0.9)	1.6487200	1.6486820	1.649067	1.66603
(1.5, 0.4, 0.3)	4.9530300	4.9530324	4.9530327	4.95316
(1.5, 0.4, 0.6)	3.6693000	3.6692940	3.6693312	3.6323
(1.5, 0.4, 0.9)	2.7182800	2.7182169	2.718851	2.74682
𝒶 = 0.9
(𝓍,𝓎,𝓉i)	$\mathcal{U}$(𝓍,𝓎,𝓉i)	LRPSM	MGMLFM [40]	FNDM [41]
(0.5, 0.4, 0.3)	1.8221188	1.7415823	1.74158	1.74178
(0.5, 0.4, 0.6)	1.3498588	1.3063222	1.3063996	1.31047
(0.5, 0.4, 0.9)	1.0000000	1.0048484	1.0058445	1.02931
(1, 0.4, 0.3)	3.0041660	2.8713838	2.871385	2.8717
(1, 0.4, 0.6)	2.2255400	2.1537611	2.1538889	2.1606
(1, 0.4, 0.9)	1.6487200	1.6567149	1.6583572	1.69705
(1.5, 0.4, 0.3)	4.9530300	4.7341116	4.7341143	4.73464
(1.5, 0.4, 0.6)	3.6693000	3.5509518	3.5511624	3.56223
(1.5, 0.4, 0.9)	2.7182800	2.7314611	2.7341689	2.79796
𝒶 = 0.75
(𝓍,𝓎,𝓉i)	$\mathcal{U}$(𝓍,𝓎,𝓉i)	LRPSM	MGMLFM [40]	FNDM [41]
(0.5, 0.4, 0.3)	1.8221188	1.6241658	1.6241897	1.6256
(0.5, 0.4, 0.6)	1.3498588	1.2598808	1.2607914	1.27791
(0.5, 0.4, 0.9)	1.0000000	1.0258447	1.033497	1.10414
(1, 0.4, 0.3)	3.0041660	2.6777967	2.6778362	2.68017
(1, 0.4, 0.6)	2.2255400	2.0771922	2.0786935	2.10692
(1, 0.4, 0.9)	1.6487200	1.6913320	1.7039486	1.82042
(1.5, 0.4, 0.3)	4.9530300	4.4149404	4.4150055	4.41885
(1.5, 0.4, 0.6)	3.6693000	3.4247111	3.4271863	3.47372
(1.5, 0.4, 0.9)	2.7182800	2.7885351	2.809336	3.00136

show a comparison of the obtained approximate solutions for the systems of time-FPDEs (7), (20) and (44), respectively with previous results generated by the existing method as MGMLFM [40], and FNDM [41] at various values of $𝓪$. As it is evident from the comparison results, the results obtained by Laplace RPSM are close to the exact solutions faster than the mentioned methods.

The 3D plots behavior of the approximate solutions of the time-FPDEs (7), (20), and (44) by Laplace RPSM are shown respectively in Figure 1, Figure 2 and Figure 3 at various values of $𝓪$ which are compared with the exact solutions on their domains. Obviously, from these figures, it can be deduced that the geometric behaviors almost agree and strongly match each other, particularly when the integer order derivative is considered. From these graphics, we can conclude that the dynamic behaviors match and correspond well with each other, specifically when the standard order derivative is considered. Moreover, Figure 4 and Figure 5 demonstrate the behavior of the obtained Laplace RPS solutions for the systems of the time-FPDEs (7) and (20) at various values of $𝓪$. It is observed from these figures that the Laplace RPSM approximate solutions match with solutions at $𝓪=1$, and this reinforces the effectiveness of the proposed method.

6. Conclusions

This investigation of time-FPDEs with initial conditions constructs a proper framework for the mathematical modeling of several fractional problems that appear in physical and engineering applications. The current work has introduced the analytical and approximate solutions for known systems of nonlinear time-FPDEs via applying Laplace RPSM. Three nonlinear time-FPDEs systems, including Broer-Kaup and Burgers’ systems, have been investigated utilizing Caputo-time fractional derivatives. The exact and the Laplace RPS solutions have been displayed numerically and graphically at various values of the fractional order 𝒶 over (0, 1]. The analysis of simulation results revealed that the Laplace RPS solutions are in imminent consistency with each other, as well as with the exact solutions at integer-order of 𝒶, which confirms the performance of the proposed method. Numerical comparisons of the obtained results with the results previously calculated by other numerical methods, such as modified generalized Mittag–Leffler function method MGMLFM [40] and fractional natural decomposition method FNDM [41], have been achieved, which indicates the high accuracy and effectiveness of the Laplace RPSM. Consequently, the analysis of attained results and their simulations confirm that the Laplace RPSM is an easy and systematic, robust, efficient, and suitable instrument to generate analytical and approximate solutions of several fractional physical and engineering problems with fewer computations and iteration steps.