1 Introduction

Nonlinear partial differential equations (PDEs) are used to model many physical phenomena arising in sciences and engineering. As a result of their considerable applications and popularity, much attention has been devoted to develop an accurate and efficient numerical method for solving PDEs. Consider a one-dimensional nonlinear parabolic partial differential equation of the form

$$ Y_{t} + \alpha {Y^{m}} {Y_{\xi }} - \beta Y_{\xi \xi }=\gamma f(Y), $$
(1)

where α, β, γ are real parameters, m is a positive integer, and \(f(Y)\) is a nonlinear function. The initial and boundary conditions are

$$\begin{aligned}& Y(\xi ,0)=Y_{0}(\xi ),\quad a\leq \xi \leq b, t> 0, \\ \end{aligned}$$
(2)
$$\begin{aligned}& Y(a,t)=g_{1}(t), \qquad Y(b,t)=g_{2}(t). \end{aligned}$$
(3)

When \(f(Y)= Y(1-Y^{m})(Y^{m}-\epsilon )\), Eq. (1) defines the generalized Burger–Huxley (GBH) equation

$$ Y_{t} + \alpha {Y^{m}} {Y_{\xi }} - \beta Y_{\xi \xi }=\gamma Y\bigl(1-Y^{m} \bigr) \bigl(Y^{m}- \epsilon \bigr),\quad 0 \leq \epsilon \leq 1, t>0. $$
(4)

Equation (4) describes the interaction between convection, reaction, and diffusion processes. When \(\alpha =0\), \(\beta =1\), and \(m=1\), Eq. (4) reduces to the Huxley equation investigating wall motion in liquid crystallography and propagation of pulse in nerve fibers. For \(\gamma =0\), \(m=1\), and \(\beta =1\), Eq. (4) reduces to the Burger equation used for analysis of nonlinear wave propagation, aspect of turbulence, traffic flows, and shock waves [1]. Similarly, when \(f(Y)=Y(1-Y^{m})\), Eq. (1) becomes the generalized Burger Fisher (GBF) equation

$$ Y_{t} + \alpha {Y^{m}} {Y_{\xi }} - \beta Y_{\xi \xi }=\gamma Y\bigl(1-Y^{m} \bigr),\quad t>0. $$
(5)

The generalized Burger–Fisher equation has a wide application in fluid mechanics, gas dynamics, plasma physics, number theory, elasticity, and heat conduction problems [2]. Equation (5) is a highly nonlinear model, which includes a combination of reaction, convection, and diffusion mechanisms. When \(\gamma =0\) and \(m=1\), Eq. (5) reduces to the Fisher equation having applications in population biology, chemistry, and biological sciences such as spreading of bacterial colonies, spread of reaction fronts in chemically bistable systems, and switching in nonlinear optics [1].

Next, we consider the following two-dimensional coupled viscous Burger equations:

$$\begin{aligned}& Y_{t}+\mu \{Y_{\xi \xi }+Y_{\eta \eta } \}+\nu (YY_{\xi })+ \alpha (YZ)_{\xi }+\gamma \{YY_{\xi }+ZY_{\eta } \}=0, \end{aligned}$$
(6)
$$\begin{aligned}& Z_{t}+\mu \{Z_{\xi \xi }+Z_{\eta \eta } \}+\nu (ZZ_{\xi })+ \beta (YZ)_{\xi }+\gamma \{YZ_{\xi }+ZZ_{\eta } \}=0, \quad (\xi , \eta )\in \Gamma , \end{aligned}$$
(7)

with initial and boundary conditions

$$\begin{aligned}& Y(\xi , \eta , t)=Y_{0}(\xi , \eta ),\qquad Z(\xi , \eta , t)=Z_{0}( \xi , \eta ), \quad (\xi , \eta )\in \Gamma , t=0, \\& Y(\xi , \eta , t)=h(t),\qquad Z(\xi , \eta , t)=g(t),\quad (\xi , \eta )\in \partial \Gamma , t\ge 0, \end{aligned}$$
(8)

where Γ and Γ represent the domain and its boundary, respectively, and μ, η, α, and β are arbitrary constants. The system was introduced by Esipov [3] to study a model of polydispersive sedimentation. This system has numerous applications in science and engineering such as gas dynamics, viscous flow of turbulence, shock waves, sedimentation of particles in fluid suspension, elasticity, and heat conduction problems [4, 5].

Many numerical methods have been applied to approximate solutions of these equations. For example, the finite difference method [6], spectral method [7], differential quadrature method [8, 9] and Adomian decomposition method [10]. Khattak [11] used a meshfree collocation method, whereas Zhu and Kang [12] applied B-spline interpolation for solution the Burger–Fisher equation. Celik [13] studied the Haar wavelet method for solving GBH equation. Dehghan [14] worked on a mixed collocation and finite difference method. Haq et al. [15] numericallly solved the Burger–Huxley equation using a meshless method of line. Zhang et al. [16] used the local discontinuous Galerkin method, whereas in [17] the author proposed a pseudospectral method for approximation of GBF equation. Wasim et al. [18] used a hybrid B-spline collocation technique for approximation of GBH and GBF equations. Mittal and Tripathe [2] proposed a cubic B-spline technique.

The modified Burger equation has been investigated in [19, 20] using a meshless method and hybrid Haar wavelet finite difference method. The author of [21] obtained approximate solution of coupled Burger equations using Adomian–Pade technique. Khater et al. [5] explained the cubic-spline collocation method for solving the coupled Burger equations. Recently, Mittal and Jiwari [22] obtained an approximate solution of one-dimensional coupled Burger equation with the help of the differential quadrature method. Dehghan et al. [23, 24] proposed a mixed finite difference and Galerkin method and multisymplectic box method for numerical study of Burgers equations. Oruc et al. [25, 26] applied a unified finite difference Chebyshev wavelet approach for time fractional Burger equations. The same authors studied the Chebyshev wavelet method for approximation of coupled Burgers equations [27]. Ali et al. [4] applied a meshfree collocation method based on the Crank–Nicolson method for time discretization and radial basis function for space discretization to solve two-dimensional coupled Burger equations, whereas the meshless method of radial basis functions (RBFs) and local RBFs were described in [28, 29] for approximate solutions of Burger-type equations. In [30] a multiscale variational algorithm was combined with the Kriging element-free Galerkin method to produce the discontinuous solutions of Burgers type equations. Srivastava et al. [31] studied a fully implicit finite difference scheme for solving two-dimensional coupled viscous Burger equations.

In this work, we compute numerical solutions of the generalized Burger–Huxley, Burger–Fisher, and coupled Burger equations using mixed Lucas and Fibonacci polynomials combined with finite differences. The main advantage of the proposed scheme is that the higher-order derivatives can be easily computed using relation of Lucas and Fibonacci polynomials. Moreover, the proposed scheme produces better accuracy for small number of collocation points, which reduces the computational cost. These polynomials have considerable applications in the area of ordinary differential equations. For example, Elhameed and Youssri [32, 33] described connection between Chebyshev and Lucas polynomials and obtained accurate solutions of boundary value problems. In [3436] the author implemented the Lucas polynomials for solutions of fractional and coupled fractional differential equations in the Caputo sense. Mostefa [37] proposed the Lucas sequence for approximation of integro-differential equations. Cetin [38] obtained numerical solution of higher-order differential equations using the Lucas polynomial approach. Farshid et al. [39] proposed a Fibonacci polynomial approach for numerical solution of Volterra–Fredholm integral differential equations. Bayku [40] presented a hybrid Taylor–Lucas polynomial method and obtained numerical solution of delay difference equations. Oruc [41, 42] for the first time applied these polynomials for solution of time-dependent partial differential equations called a mixed Lucas and Fibonacci polynomial technique.

The rest of the paper is organized as follows. In Sect. 2, we discuss description of solution methodology. In Sect. 3, describe the stability of the method. In Sects. 4 and 5, we present numerical experiments followed by conclusion of the paper.

2 Methodology

In this section, we give a description of the proposed method for two different cases. The suggested technique will be tested by some examples.

Case.1 Nonlinear PDEs

For solution of Eq. (1) using the proposed technique, we discretize the time derivative of the equation with finite differences and apply the θ-weighted scheme to its spatial part to get

$$\begin{aligned} &\frac{1}{\delta t}{\bigl[Y^{n+1}-Y^{n} \bigr]}+\theta \bigl\{ \alpha \bigl({Y^{m}} {Y_{\xi }} \bigr)^{n+1} - \beta Y_{\xi \xi }^{n+1}-\gamma \bigl( f(Y)\bigr)^{n+1} \bigr\} \\ &\quad {}+(1-\theta ) \bigl\{ \alpha {Y^{m}} {Y_{\xi }} - \beta Y_{\xi \xi }- \gamma f(Y) \bigr\} ^{n}=0,\quad 0\leq \theta \leq 1, \end{aligned}$$
(9)

where \(Y^{n}=Y(\xi ,t^{n})\), \(t^{n}=n\delta t\), \(n=1,2,\dots ,N\), δt is the time step. The nonlinear term in Eq. (9) is linearized using the lagging method given as

$$ \bigl(Y^{m}Y_{\xi }\bigr)^{n+1}= \bigl(Y^{n}\bigr)^{m}Y_{\xi }^{n+1}. $$
(10)

Using Eq. (10) in Eq. (9), we get

$$\begin{aligned}& Y^{n+1}+\theta \delta t \bigl\{ \alpha \bigl(Y^{n}\bigr)^{m}Y_{\xi }^{n+1} - \beta Y_{\xi \xi }^{n+1}-\gamma \bigl( f(Y)\bigr)^{n+1} \bigr\} \\& \quad =Y^{n}+( \theta -1) \delta t \bigl\{ \alpha \bigl({Y^{m}} {Y_{\xi }}\bigr)^{n} - \beta Y_{ \xi \xi }^{n}-\gamma \bigl(f(Y)\bigr)^{n} \bigr\} . \end{aligned}$$
(11)

Approximating \(Y^{n}(\xi )\) by Lucas polynomials is as follows:

$$ {Y}^{n}{(\xi )}=\sum_{k=1}^{N}C_{k}^{n}L_{k}( \xi )=WC^{n}, $$
(12)

where \(C_{k}^{n}\) are unknown coefficients to be computed, and \(L_{k}(\xi )\) are the Lucas polynomials defined by [41]

$$ L_{k}(\xi )=\xi L_{k-1}(\xi )+L_{k-2}(\xi ),\quad k\geq 2, \quad \text{with}\quad L_{o}(\xi )=2,\qquad L_{1}(\xi )= \xi . $$

To determine \(C_{k}^{n}\), we use the collocation method. At collocation points \(\xi _{i}=a+i \delta \xi \), Eq. (12) can be written as

$$ Y^{n}(\xi _{i})=\sum _{k=0}^{N}C_{k}^{n}L_{k}( \xi _{i}),\quad i=1, \dots ,N. $$
(13)

Putting the values from Eq. (13) into Eq. (11), we obtain the following system of N equations:

$$\begin{aligned}& {\sum_{k=1}^{N} C^{n+1}_{k} L_{k}(\xi _{i})} \\& \qquad {}+{\theta } {\delta t} { \Biggl[\alpha \Biggl\{ \sum_{k=1}^{N} C^{n}_{k} L_{k}(\xi _{i}) \Biggr\} ^{m}{\sum_{k=1}^{N}C_{k}^{n+1}L_{k}^{\prime }( \xi _{i})}-\beta { \sum_{k=1}^{N}C_{k}^{n+1}L_{k}^{\prime \prime }( \xi _{i})}-\gamma \bigl( f(Y)\bigr)^{n+1} \Biggr]} \\& \quad ={\sum_{k=1}^{N} C^{n}_{k} L_{k}(\xi _{i})} \\& \qquad {}+{( \theta -1)} {\delta t} { \Biggl[\alpha \Biggl\{ \sum_{k=1}^{N} C^{n}_{k} L_{k}(\xi _{i}) \Biggr\} ^{m}{\sum_{k=1}^{N}C_{k}^{n}L_{k}^{\prime }( \xi _{i})}-\beta { \sum_{k=1}^{N}C_{k}^{n}L_{k}^{\prime \prime }( \xi _{i})}-\gamma \bigl( f(Y)\bigr)^{n} \Biggr]}. \end{aligned}$$
(14)

The primes in Eq. (14) represent differentiation with respect to ξ, which allows us to replace \(L_{k}\) by Fibonacci polynomials [41]:

$$ L_{k}^{\prime }(\xi )=kF_{k}( \xi ),\qquad L_{k}^{\prime \prime }(\xi )=kF_{k}(\xi )D, $$
(15)

where \(F_{k}(\xi )\) are the Fibonacci polynomials defined as [41].

$$ F_{k}(\xi )=\xi F_{k-1}(\xi )+F_{k-2}(\xi ) \quad \text{for } k \geq {2}, \quad \text{with}\quad F_{0}(\xi )=0,\qquad F_{1}(\xi )=1, $$

and D is differentiation matrix given by [41]

D= [ 0 0 0 0 d 0 ] ,

where d is the square matrix of order N defined by

$$ d_{m,n}= \textstyle\begin{cases} m(-1)^{\frac{(n-m-1)}{2}} & \text{if }n>m, n-m\mbox{ odd}, \\ 0 & \text{otherwise}. \end{cases} $$

Substituting the values from Eq. (15) into Eq. (14), we can write

$$\begin{aligned}& {\sum_{k=1}^{N} C^{n+1}_{k} L_{k}(\xi _{i})} \\& \qquad {}+{ \theta } {\delta t} { \Biggl[\alpha \Biggl\{ \sum_{k=1}^{N} C^{n}_{k} L_{k}(\xi _{i}) \Biggr\} ^{m}{\sum_{k=1}^{N}C_{k}^{n+1}kF_{k}( \xi _{i})}-\beta { \sum_{k=1}^{N}C_{k}^{n+1}kF_{k}( \xi _{i})D}-\gamma \bigl( f(Y)\bigr)^{n+1} \Biggr]} \\& \quad = {\sum_{k=1}^{N} C^{n}_{k} L_{k}(\xi _{i})}+{( \theta -1)} {\delta t} \Biggl[\alpha \Biggl\{ \sum_{k=1}^{N} C^{n}_{k} L_{k}(\xi _{i}) \Biggr\} ^{m}{\sum_{k=1}^{N}C_{k}^{n}kF_{k}( \xi _{i})} \\& \qquad {}-\beta {\sum_{k=1}^{N}C_{k}^{n}kF_{k}( \xi _{i})D}-\gamma \bigl( f(Y)\bigr)^{n} \Biggr]. \end{aligned}$$
(16)

Similarly, Eqs. (2) and (3) can be transformed to

$$ \begin{aligned} &\sum_{k=1}^{N} C^{1}_{k} L_{k}(\xi _{i}))=Y_{0}(\xi _{i}), \qquad \sum _{k=1}^{N}C_{k}^{n+1}L_{k}( \xi _{1})=g_{1}(t), \\ &\sum_{k=1}^{N}C_{k}^{n+1}L_{k}( \xi _{N})=g_{2}(t)\quad i=1,\dots ,N. \end{aligned}$$
(17)

In matrix form, Eqs. (16)–(17) can be written as

$$ HC^{n+1}=GC^{n}+A^{n+1}, $$
(18)

where H, G, and A are square matrices of order N with components given by

$$\begin{aligned}& H_{ik}= \textstyle\begin{cases} L_{k}(\xi _{i})+\delta t \theta \{ \alpha L_{k}(\xi _{i}) kF_{k}( \xi _{i}) -\beta kF_{k}(\xi _{i})D-\gamma f(Y) \} , \\ \hphantom{L_{k}(\xi _{i}),}\quad i={2, \dots ,N-1}, k=1,2,\ldots,N, \\ L_{k}(\xi _{i}), \quad i={1,N}, k=1,2,\ldots,N, \end{cases}\displaystyle \end{aligned}$$
(19)
$$\begin{aligned}& G_{ik}= \textstyle\begin{cases} L_{k}(\xi _{i})+\delta t (\theta -1) \{ \alpha L_{k}(\xi _{i}) kF_{k}( \xi _{i}) -\beta kF_{k}(\xi _{i})D-\gamma f(Y) \} , \\ \hphantom{0,} \quad i=2, \dots ,N-1 , k=1,2,\dots ,N, \\ 0, \quad i=1,N, k=1,2,\dots ,N, \end{cases}\displaystyle \end{aligned}$$
(20)
$$\begin{aligned}& A_{ik}= \textstyle\begin{cases} 0, &i=2,\ldots,N-1 , k=1,2,\ldots,N, \\ g^{n+1}(\xi _{i}), &i=1,N , k=1,2,\ldots,N. \end{cases}\displaystyle \end{aligned}$$
(21)

Solution of Eq. (18) give required unknowns C, and hence a solution of problem (1) can be obtained with the help of Eq. (13).

Case 2: Coupled PDEs

To construct a scheme for coupled Burger Eqs. (6)–(7), discretizing the temporal and spatial parts in a similar way as discussed in case 1, we have

$$\begin{aligned} &\frac{Y^{n+1}-Y^{n}}{\delta t}+\theta \bigl[\mu \bigl(Y_{\xi \xi }^{n+1}+Y_{ \eta \eta }^{n+1} \bigr)+\nu (YY_{\xi })^{n+1}+\alpha \bigl((YZ_{\xi })^{n+1}+(ZY_{\xi })^{n+1} \bigr) \\ &\qquad {}+\gamma { \bigl[(YY_{\xi })^{n+1}+(ZY_{\eta })^{n+1} \bigr]} \bigr] \\ &\quad =(\theta -1) \bigl[\mu \bigl(Y_{\xi \xi }^{n} +Y_{\eta \eta }^{n}\bigr)+\nu (YY_{\xi })^{n}+ \alpha \bigl((YZ_{\xi })^{n}+(ZY_{\xi })^{n} \bigr) \\ &\qquad {}+\gamma \bigl[(YY_{\xi })^{n}+(ZY_{ \eta })^{n} \bigr] \bigr], \end{aligned}$$
(22)
$$\begin{aligned} &\frac{Z^{n+1}-Z^{n}}{\delta t} +\theta \bigl[\mu \bigl(Z_{\xi \xi }^{n+1}+Z_{ \eta \eta }^{n+1} \bigr)+\nu (ZZ_{\xi })^{n+1}+\beta \bigl((YZ_{\xi })^{n+1}+(ZY_{\xi })^{n+1} \bigr) \\ &\qquad {}+ \gamma { \bigl[(YZ_{\xi })^{n+1}+(ZZ_{\eta })^{n+1} \bigr]} \bigr] \\ &\quad =(\theta -1) \bigl[\mu \bigl(Z_{\xi \xi }^{n}+Z_{\eta \eta }^{n} \bigr)+\nu (ZZ_{\xi })^{n}+\beta \bigl((YZ_{\xi })^{n}+(ZY_{\xi })^{n} \bigr) \\ &\qquad {}+\gamma { \bigl[(YZ_{\xi })^{n}+(ZZ_{ \eta })^{n} \bigr]} \bigr]. \end{aligned}$$
(23)

For \(\theta =1/2\), these equations become the well-known Cran–Nicolson scheme with accuracy \(O(\delta t^{2})\) [43]. For nonlinear terms, we use the lagging method

$$ (YY_{\xi })^{n+1}=Y^{n+1}Y_{\xi }^{n}. $$
(24)

Using Eq. (24) in Eqs. (22)–(23) and waiving the error terms, we get

$$\begin{aligned} &Y^{n+1}+\delta t \theta \bigl[\mu \bigl(Y_{\xi \xi }^{n+1}+Y_{\eta \eta }^{n+1} \bigr)+ \nu Y^{n+1}Y_{\xi }^{n}+\alpha \bigl(Y^{n+1}Z_{\xi }^{n}+Z^{n}Y_{\xi }^{n+1} \bigr) \\ &\qquad {}+ \gamma \bigl(Y^{n+1}Y_{\xi }^{n}+Z^{n}Y_{\eta }^{n+1} \bigr) \bigr] \\ &\quad =Y^{n}+(\theta -1) \delta t \bigl[\mu \bigl(Y_{\xi \xi }^{n}+Y_{\eta \eta }^{n} \bigr)+ \nu (YY_{\xi })^{n}+\alpha \bigl((YZ_{\xi })^{n}+(ZY_{\xi })^{n} \bigr) \\ &\qquad {}+\gamma (YY_{ \xi }+ZY_{\eta })^{n} \bigr], \end{aligned}$$
(25)
$$\begin{aligned} &Z^{n+1}+\delta t \theta \bigl[\mu \bigl(Z_{\xi \xi }^{n+1}+Z_{\eta \eta }^{n+1} \bigr)+ \nu Z^{n+1}Z_{\xi }^{n}+\beta \bigl(Y^{n}Z_{\xi }^{n+1}+Z^{n+1}Y_{\xi }^{n} \bigr) \\ &\qquad {}+ \gamma \bigl(Y^{n}Z_{\xi }^{n+1}+Z^{n+1}Z_{\eta }^{n} \bigr) \bigr] \\ &\quad = Z^{n}+(\theta -1) \delta t \bigl[\mu \bigl(Z_{\xi \xi }^{n}+Z_{\eta \eta }^{n} \bigr)+\nu (ZZ_{\xi })^{n}+\beta \bigl((YZ_{\xi })^{n}+(ZY_{\xi })^{n} \bigr) \\ &\qquad {}+ \gamma \bigl((YZ_{\xi })^{n}+(ZZ_{\eta })^{n} \bigr) \bigr]. \end{aligned}$$
(26)

Now we approximate \(Y^{n}\) and \(Z^{n}\) by Lucas polynomials as follows:

$$ Y^{n}(\xi , \eta )=\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j}),\qquad Z^{n}(\xi , \eta )=\sum_{k=1}^{N} \sum_{m=1}^{N} \lambda _{km}^{n}L_{k}(\xi _{i})L_{m}( \eta _{j}), $$
(27)

where \(\lambda _{km}^{n}\) and \(C_{km}^{n}\) are unknown coefficients, and \({\xi _{i}=\eta _{i}=a+(i-1)\,d\xi }_{i}\), \(d{\xi }=d{\eta }\), are the regular collocation points or the Chebyshev–Gauss–Lobatto (CGL) collocation points

$$ {\xi _{i}=\eta _{i}=a+\frac{b-a}{2}\bigl(1- \cos \bigl((i-1)\pi /M\bigr)\bigr)} $$

with \(a=\xi _{1}=\eta _{1}\) and \(b=\xi _{M}=\eta _{M}\); is spacial step size. Plugging Eq. (27) into Eqs. (25)–(26), we get

$$\begin{aligned}& \sum_{k=1}^{N} \sum_{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}( \eta _{j}) \\& \qquad {}+{\delta t} { \theta } \Biggl\{ \mu \Biggl( \sum_{k=1}^{N} \sum_{m=1}^{N}C_{km}^{n+1}L_{k}^{\prime \prime }( \xi _{i})L_{m}({\eta _{j}})+\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}^{\prime \prime }({\eta _{j}}) \Biggr) \\& \qquad {} +{\nu }\sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j}) \\& \qquad {}+\alpha \sum _{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}(\eta _{j}) \\& \qquad {}* \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j})+ \alpha \sum _{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}( \xi _{i})L_{m}(\eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n+1}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j}) \\& \qquad {}+ \gamma \sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j}) \\& \qquad {} +\gamma \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}^{\prime }( \eta _{j}) \Biggr\} \\& \quad =\sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}(\eta _{j}) \end{aligned}$$
(28)
$$\begin{aligned}& \qquad {} +{\delta t} {(\theta -1)} \Biggl\{ \mu \Biggl( \sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}^{\prime \prime }( \xi _{i})L_{m}({\eta _{j}})+\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}^{\prime \prime }({\eta _{j}}) \Biggr) \\& \qquad {} +{\nu }\sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j})+\alpha \sum _{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j}) \\& \qquad {}* \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j})+ \alpha \sum _{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}( \xi _{i})L_{m}(\eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j}) \\& \qquad {}+\gamma \sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j}) \\& \qquad {} + \gamma \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}^{\prime }( \eta _{j}) \Biggr\} , \\& \sum_{k=1}^{N} \sum_{m=1}^{N}\lambda _{km}^{n+1}L_{k}(\xi _{i})L_{m}( \eta _{j}) \\& \qquad {}+{\delta t} {\theta } \Biggl\{ \mu \Biggl( \sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n+1}L_{k}^{\prime \prime }(\xi _{i})L_{m}({\eta _{j}})+\sum _{k=1}^{N} \sum_{m=1}^{N} \lambda _{km}^{n+1}L_{k}(\xi _{i})L_{m}^{\prime \prime }({\eta _{j}}) \Biggr) \\& \qquad {} +{\nu }\sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n+1}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j})+\beta \sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n+1}L_{k}(\xi _{i})L_{m}(\eta _{j}) \\& \qquad {}* \sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j})+ \beta \sum _{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n+1}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j}) \\& \qquad {} +\gamma \sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n+1}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j}) \\& \qquad {}+\gamma \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n+1}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}(\xi _{i})L_{m}^{\prime }( \eta _{j}) \Biggr\} \\& \quad = \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})+{\delta t} {( \theta -1)} \\& \qquad {}* \Biggl\{ \mu \Biggl( \sum_{k=1}^{N} \sum_{m=1}^{N}\lambda _{km}^{n}L_{k}^{\prime \prime }(\xi _{i})L_{m}({\eta _{j}})+ \sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}(\xi _{i})L_{m}^{\prime \prime }({ \eta _{j}}) \Biggr) \\& \qquad {} +{\nu }\sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j})+\beta \sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}(\xi _{i})L_{m}(\eta _{j}) \\& \qquad {}* \sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}( \eta _{j})+ \beta \sum _{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j}) \\& \qquad {}+\gamma \sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}^{\prime }( \xi _{i})L_{m}(\eta _{j}) \\& \qquad {}+\gamma \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n}L_{k}( \xi _{i})L_{m}( \eta _{j})\sum _{k=1}^{N}\sum_{m=1}^{N} \lambda _{km}^{n}L_{k}(\xi _{i})L_{m}^{\prime }( \eta _{j}) \Biggr\} . \end{aligned}$$
(29)

In these equations, “*” represents the usual product. Similarly, boundary conditions given in Eq. (8) take the form

$$\begin{aligned}& \sum_{k=1}^{N}\sum _{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}(\eta _{j})=h_{1} \bigl(t^{n+1}\bigr),\qquad \sum_{k=1}^{N} \sum_{m=1}^{N}C_{km}^{n+1}L_{k}( \xi _{i})L_{m}( \eta _{j})=h_{2} \bigl(t^{n+1}\bigr), \end{aligned}$$
(30)
$$\begin{aligned}& \sum_{k=1}^{N}\sum _{m=1}^{N}\lambda _{km}^{n+1}L_{k}( \xi _{i})L_{m}( \eta _{j})=g_{1} \bigl(t^{n+1}\bigr),\qquad \sum_{k=1}^{N} \sum_{m=1}^{N} \lambda _{km}^{n+1}L_{k}(\xi _{i})L_{m}( \eta _{j})=g_{2}\bigl(t^{n+1}\bigr), \end{aligned}$$
(31)

considering the relation between the Lucas and Fibonacci polynomials

$$ L_{k}^{\prime }(\xi )=kF_{k}( \xi ),\qquad L_{k}^{\prime \prime }(\xi )=kF_{k}(\xi )D, $$
(32)

where \(F(\xi )\) and D have the same meaning as before. Putting values from Eq. (32) into Eqs. (28)–(29), the matrix form of Eqs. (28)–(31) can be written as

$$\begin{aligned}& A+{\delta t} {\theta } \bigl\{ \mu (Q_{3}+Q_{4})+\nu A*Y_{\xi }^{n}+ \alpha \bigl(A*Z_{\xi }^{n}+Q_{1}*Z^{n} \bigr) \end{aligned}$$
(33)
$$\begin{aligned}& \qquad {}+\gamma \bigl(A*Y_{ \xi }^{n}+Q_{2}*Z^{n} \bigr) \bigr\} C^{n+1} \\& \quad =A+{\delta t} {(\theta -1)} \bigl\{ \mu (Q_{3}+Q_{4})+ \nu A*Y_{\xi }^{n}+ \alpha \bigl(A*Z_{\xi }^{n}+Q_{1}*Z^{n} \bigr) \\& \qquad {}+\gamma \bigl(A*Y_{ \xi }^{n}+Q_{2}*Z^{n} \bigr) \bigr\} C^{n}+H^{n+1}, \\& A+{\delta t} {\theta } \bigl\{ \mu (Q_{3}+Q_{4})+\nu A*Z_{\xi }^{n}+ \beta \bigl(A*Y_{\xi }^{n}+Q_{1}*Y^{n} \bigr) \\& \qquad {}+\gamma \bigl(Q_{1}*Y^{n}+A*Z_{ \eta }^{n} \bigr) \bigr\} \lambda ^{n+1} \\& \quad =A+{\delta t} {(\theta -1)} \bigl\{ \mu (Q_{3}+Q_{4})+ \nu A*Z_{\xi }^{n}+ \beta \bigl(A*Y_{\xi }^{n}+Q_{1}*Y^{n} \bigr) \\& \qquad {}+\gamma \bigl(Q_{1}*Y^{n}+A*Z_{ \eta }^{n} \bigr) \bigr\} \lambda ^{n}+G^{n+1}, \end{aligned}$$
(34)

where

$$\begin{aligned}& A=L_{k}(\xi _{i})L_{m}(\eta _{j})\}_{k,m=1}^{N},\qquad Q_{1}=kF_{k}( \xi _{i})L_{m}(\eta _{j}) \}_{k,m=1}^{N},\qquad Q_{2}=L_{k}(\xi _{i})mF_{m}( \eta _{j}) \}_{k,m=1}^{N}, \\& Q_{3}=kF_{k}(\xi _{i})DL_{m}( \eta _{j})\}_{k,m=1}^{N},\qquad Q_{4}=L_{k}( \xi _{i})mF_{m}D(\eta _{j}) \}_{k,m=1}^{N}, \end{aligned}$$

and

$$\begin{aligned}& H^{n+1}= \bigl\{ h_{1}^{n+1},0,0, \ldots,h_{2}^{n+1} \bigr\} ^{T}, \\& G^{n+1}= \bigl\{ g_{1}^{n+1},0,0, \ldots,g_{2}^{n+1} \bigr\} ^{T}. \end{aligned}$$

Equations (33)–(34) can be written as

$$\begin{aligned}& C^{n+1}=M_{1}^{-1}N_{1}C_{1}^{n}+M_{1}^{-1}H^{n+1}, \end{aligned}$$
(35)
$$\begin{aligned}& \lambda ^{n+1}=M_{2}^{-1}N_{2}C_{2}^{n}+M_{2}^{-1}G^{n+1}, \end{aligned}$$
(36)

where

$$\begin{aligned}& M_{1}=A+{\delta t} {\theta } { \bigl\{ \mu (Q_{3}+Q_{4})+ \nu A*Y_{\xi }^{n}+\alpha \bigl(A*Z_{\xi }^{n}+Q_{1}*Z^{n} \bigr)+ \gamma \bigl(A*Y_{\xi }^{n}+Q_{2}*Z^{n} \bigr) \bigr\} }, \\& N_{1}=A+{\delta t} {(\theta -1)} \bigl\{ \mu (Q_{3}+Q_{4})+\nu A*Y_{ \xi }^{n}+ \alpha \bigl(A*Z_{\xi }^{n}+Q_{1}*Z^{n} \bigr) \\& \hphantom{N_{1}=}{}+\gamma \bigl(A*Y_{\xi }^{n}+Q_{2}*Z^{n} \bigr) \bigr\} , \\& M_{2}=A+{\delta t} {\theta } { \bigl\{ \mu (Q_{3}+Q_{4})+ \nu A*Z_{\xi }^{n}+\beta \bigl(A*Y_{\xi }^{n}+Q_{1}*Y^{n} \bigr)+ \gamma \bigl(Q_{1}*Y^{n}+A*Z_{\eta }^{n} \bigr) \bigr\} }, \\& N_{2}=A+{\delta t} {(\theta -1)} \bigl\{ \mu (Q_{3}+Q_{4})+\nu A*Z_{ \xi }^{n}+ \beta \bigl(A*Y_{\xi }^{n}+Q_{1}*Y^{n} \bigr) \\& \hphantom{N_{2}=}{}+\gamma \bigl(Q_{1}*Y^{n}+A*Z_{\eta }^{n} \bigr) \bigr\} . \end{aligned}$$

Since \(Y^{n}=AC^{n}\) and \(Z^{n}=A\lambda ^{n}\), we get

$$\begin{aligned}& Y^{n+1}=AM_{1}^{-1}N_{1}A^{-1}Y^{n}+AM_{1}^{-1}H^{n+1}, \end{aligned}$$
(37)
$$\begin{aligned}& Z^{n+1}=AM_{2}^{-1}N_{2}A^{-1}Z^{n}+AM_{2}^{-1}G^{n+1}. \end{aligned}$$
(38)

If \(M_{1}\), \(M_{2}\) are fully ranked, then \(M_{1}^{-1}\) and \(M_{2}^{-1}\) exist [44, 45]. In our computation, these matrices are fully ranked, which is checked using Matlab command \(\operatorname{rank}(M_{i})\). Therefore system (35)–(36) can be solved for unknown \(C^{n}\) and \(\lambda ^{n}\), and a solution of original coupled equations can be obtained from Eq. (27).

3 Numerical stability analysis

To check the stability of the proposed technique, we use the matrix method. For this purpose, first rewrite Eq. (18) as follows:

$$ C^{n+1}=H^{-1}GC^{n}+H^{-1}A^{n+1} \quad \text{for } n\geq 0. $$
(39)

If Y denotes the approximate solution, and u denotes the exact one, then the error is defined as

$$ E^{n}=u^{n+1}-Y^{n+1} \quad \text{with } Y^{n+1}=WC^{n+1}. $$
(40)

Substituting the values from Eq. (39) into Eq. (40), we get

$$ E^{n+1}=WH^{-1}GW^{-1}E^{n}=ME^{n}. $$
(41)

Here \(\textbf{M}=WH^{-1}GW^{-1}\) is known as the implication matrix. Scheme (39) is stable if the matrix M satisfies the Lax–Richtmyer stability condition

$$ \Vert {\textbf{M}} \Vert \leq 1. $$
(42)

If \(\rho(\textbf{M})\) denotes the spectral radius of a matrix M, then \(\rho (\textbf{M})\leq \|\textbf{M} \|\). The validation of this condition is described with the help of different numerical examples in the next section.

4 Numerical examples

In this section, we apply the proposed technique to some problems. Efficiency of the scheme is tested by comparing the obtained results with exact and approximate solutions available in the literature. We study the accuracy of the method in the form of the root mean square (RMS), \(L_{2}\) and \(L_{\infty }\) norms, and the time convergence rate given by

$$\begin{aligned}& L_{\infty }=\max \bigl\vert E(\xi _{i},t)-Y(\xi _{i},t) \bigr\vert _{i=1}^{N}, \qquad RMS=\sqrt{ \frac{\sum_{i=1}^{N} \vert E(\xi _{i},t)-Y(\xi _{i},t) \vert ^{2}}{N}}, \\& L_{2}=\sqrt{ \frac{ \sum_{i=1}^{N} \vert E(\xi _{i},t)-Y(\xi _{i},t) \vert ^{2}}{d\xi }},\qquad Rate= \frac{\log (\frac{L_{\infty }(2\delta t,d\xi )}{L_{\infty }( \delta t, d \xi )} )}{\log (\frac{2\delta t}{\delta t} )}. \end{aligned}$$

Example 1

Consider Eq. (4) with \(\beta =1\) and exact solution [2]

$$ Y(\xi ,t)= \biggl[\frac{\epsilon }{2}+\frac{\epsilon }{2} \tanh \bigl(a_{1}( \xi -a_{2}t)\bigr) \biggr]^{1/m},\quad t\geq 0, $$
(43)

where

$$\begin{aligned}& a_{1}=\frac{-\alpha m + m\sqrt{\alpha ^{2}+4 \gamma (1+m)}}{4(1+m)} \epsilon , \\& a_{2}=\frac{\alpha \epsilon }{(1+m)}- \frac{(1+m-\epsilon )(-\alpha +\sqrt{\alpha ^{2}+4 \gamma (1+m)})}{2(1+m)}. \end{aligned}$$

We obtain the approximate solution the proposed method taking special domains \([0, 1]\) and \([-10, 20]\). Initial and boundary conditions are taken from the exact solution. In this example, we discuss various cases for different values of the parameters α, γ, ϵ, and m appearing in Eq. (4). We compute the solution for both regular and CGL collocation points. We compare the computed results in the form of error norms with those available in the literature [2, 7]. From comparison it is clear that the present method gives better accuracy or comparable results with those available in the literature. We can observe from the comparison that the proposed method produces slightly more accurate results for CGL collocations points than for regular points. We also report and show in tables the comparison carried out for stability of the scheme in the form of spectral radius.

Case 1.1: In this case, we take \(\alpha =\gamma =1\), \(m=2\), and \(\epsilon =0.5\) with nodal points \(N=10\). The solution is computed over the domain \([0, 1]\) for different values of \(T=15, 30, 60\), and 120 with step size 0.01 and is shown in Tables 1 and 3. From Table 1 it is clear that the accuracy of the solution increases with time where the system remains stable, that is, \(\rho (\textbf{M}) < 1\). The computed results are compared with those of the collocation cubic B-spline method [2]. From the comparison it is obvious that the present technique gives better accuracy than the collocation cubic B-spline method.

Table 1 Error norms and spectral radius for the solution of Example 1 in Case 1.1

Case 1.2: In this case, we take the parameters are \(\alpha =\gamma =m=1\) and \(\epsilon =0.5\) with nodal points \(N=10\). The solution has been computed at different time levels \(T=15, 30, 60\), and 120 with step size \(\delta t=0.01\). For accuracy of the scheme, different error norms were calculated and compared with the error norms of cubic B-spline [2] and tabulated in Tables 2 and 3. From Table 2 we observe that for small times, the results of cubic B-spline are better than those of the present method, but as time increases, the accuracy of the proposed technique increases. We can also observe from the table that the spectral radius \(\rho (\textbf{M})<1\) for all time levels, which shows the stability of the scheme.

Table 2 Error norms and spectral radius for the solution of Example 1 in Case 1.2
Table 3 Error norms of Example 1 using regular points

Case 1.3: In this case, the solution was computed using CGL grid points for \(\alpha =\gamma =1\) and \(\epsilon =0.001\), and various values of \(m=1, 4, 8\). The spectral radius and error norms were computed at different values of \(T=0.2, 1\) and compared with the results available in the literature, which are shown in Table 4. From the table it is clear that the present method gives an excellent solution in comparison to available techniques. We can observe from the table that the present method gives the same accuracy for all values of the parameter m, whereas the accuracy of spectral and cubic B-spline methods [2, 7] suddenly decreases, which reflects the feasibility of the proposed scheme. The solution is also been computed using regular grid point and shown in Table 6.

Table 4 Error norms and spectral radius for the solution of Example 1 in Case 1.3

Case 1.4: In this case, we take \(\alpha =\gamma =m=1\) and \(\epsilon =0.01\). The solution is computed over the domain \([-10, 20]\) with step size \(\delta t=0.003\). The error norms and spectral radius are computed for different time levels \(T=0.1, 1, 5\), and 10, which are presented in Tables 5 and 6. From the table we can observe that the present method gives better results than those available in the literature even for large space and time domains. We can easily understand from these tables that as the domain increases, the spectral radius remains \(\rho (\textbf{M})<1\), which shows the stability of the proposed scheme. Overall, it is obvious that the present method gives better results and is flexible to implement.

Table 5 Error norms and spectral radius for the solution of Example 1 in Case 1.4
Table 6 Error norms of Example 1 using regular points

Example 2

Consider Eq. (5) with \(\beta =1\). In this case, the exact solution [2] is given by

$$ Y(\xi ,t)= \biggl[\frac{1}{2}+\frac{1}{2} \tanh \bigl(w_{1}(\xi -w_{2}t)\bigr) \biggr]^{1/m},\quad t\geq 0, $$
(44)

where

$$ w_{1}=\frac{-\alpha m}{2(1+m)}, \qquad w_{2}=\frac{\alpha }{(1+m)}+ \frac{\gamma (1+m)}{\alpha }. $$

The initial and boundary conditions are taken from the exact solution, and numerical solution is computed using the suggested technique with domain \([-1, 1]\). The error norms were computed for comparison of the proposed scheme with exact and available solutions in the literature. Two different cases were discussed for different values of the real parameter m, whereas \(\alpha =\gamma =\beta =1\).

Case 2.1: In this case, the parameter \(m=1\) with step size \(\delta t=0.0001\). The solution was obtained for \(T= 0.5, 1, 2\), and 4. The computed results are compared with the exact and available numerical solution in the form of error norms. The stability of the scheme was studied in the form of spectral radius shown in Tables 7 and 9. From Table 7 we can notice that the proposed scheme is stable in the given domain and gives a better accuracy than the collocation cubic B-spline method [2].

Table 7 Error norms and spectral radius for the solution of Example 2 in Case 2.1

Case 2.2: Now we take \(m=2\) and compute the solution for different time levels \(T=5, 10, 15\), and 20. The error norms along with spectral radius were obtained at each time level, and the results were compared with existing numerical solutions described in Tables 8 and 9. It is obvious from the table that the accuracy increases with time and converges toward the true solution. It is also clear from the table that in all the three error norms the proposed method gives excellent results at each time level as compared to available results in the literature. The table also shows that the value of spectral radius remains less than 1, which clarifies the stability of the scheme.

Table 8 Error norms and spectral radius for the solution of Example 2 in Case 2.2
Table 9 Error norms of Example 2 using regular points

Example 3

Putting \(f(Y)=Y(1-Y)(Y-0.5)\) and \(\alpha =\gamma =m=1\) in Eq. (1), we have

$$ Y_{t}+YY_{\xi }-\beta Y_{\xi \xi }=Y(1-Y) (Y-0.5) $$
(45)

with initial and boundary conditions [2]

$$\begin{aligned}& Y(\xi ,0)=\sin (\pi \xi ), \quad 0 \leq x \leq 1, \end{aligned}$$
(46)
$$\begin{aligned}& Y(0,t)=Y(1,t)=0, \quad t \ge 0. \end{aligned}$$
(47)

The solution is computed for different time levels \(T=0.1, 0.3, 0.6\), and 0.9 and various values of \(\beta = 2^{-1}, 2^{-4}, 2^{-6}\), and 2−9, respectively. The solution profile for different time levels is plotted in Fig. 1, which shows the same pattern reported in [2]. Due to unavailability of exact solution, we studied the accuracy via the stability of the method and noticed that as the value of β approaches zero, the solution diverges, and the scheme becomes unstable. From Figs. 1(a)–1(c) it is clear that the solution profile shows a proper behavior and the scheme is stable, that is, \(\rho (\textbf{M})<1\), and as the spectral radius increases, the scheme goes to an unstable region, and the solution diverges, as shown in Fig. 1(d). We can also observe from the figure that for a fixed value of β, the graph tends to zero as time increases, and for various values of β, the curve propagates to the right and follows a sharp decay. Thus the present method shows the proper behavior of Eq. (45) for various values of β and t. The same behavior of this equation was illustrated by Mahanty and Sharma [1] and Mittal [2].

Figure 1
figure 1

Approximate solutions of Example 3: (a) solution of \(\beta =2^{-1}\) with \(\rho (\textbf{M})<1\), (b) solution of \(\beta =2^{-4}\) with \(\rho (\textbf{M})<1\), (c) solution of \(\beta =2^{-6}\) with \(\rho (\textbf{M})<1\), and (d) solution of \(\beta =2^{-9}\) with \(\rho (\textbf{M})>1\)

Example 4

Consider Eq. (1) with \(\gamma =0\) and \(\alpha =m=1\) to obtain

$$ Y_{t}+YY_{\xi }-\beta Y_{\xi \xi }=0 $$
(48)

with initial and boundary conditions [1]

$$\begin{aligned}& Y(\xi ,0)=\xi \bigl(1-{\xi }^{2}\bigr),\quad 0 \leq x \leq 1, \end{aligned}$$
(49)
$$\begin{aligned}& Y(0,t)=Y(1,t)=0,\quad t \ge 0. \end{aligned}$$
(50)

We computed the solution of the problem for various values of \(\beta =2^{-3}, 2^{-6}, 2^{-9}, 2^{-11}\) and time \(T=0.1,0.3,0.6,0.9\) with step size \(\delta t=0.001\). The obtained results are presented in Fig. 2, which clarifies the behavior of the problem. Similarly to Example 3, We can notice that as the value of β decreases, the solution diverges, and the scheme become unstable. The same pattern of the problem has been reported by [8] and [1]. The accuracy of solution was discussed by means of the stability of the scheme.

Figure 2
figure 2

Approximate solutions of Example 4: (a) solution of \(\beta =2^{-3}\) with \(\rho (\textbf{M})<1\), (b) solution of \(\beta =2^{-6}\) with \(\rho (\textbf{M})<1\), (c) solution of \(\beta =2^{-9}\) with \(\rho (\textbf{M})=1.003\), and (d) solution of \(\beta =2^{-11}\) with \(\rho (\textbf{M})=1.75\)

Example 5

In this case, we consider the coupled one-dimensional Burger Eq. (7) by taking \(\mu =-1\), \(\nu =2\), and \(\gamma =0\), which leads to

$$\begin{aligned}& Y_{t}-Y_{\xi \xi }+2YY_{\xi }+ \alpha (YZ)_{\xi }=0, \\ \end{aligned}$$
(51)
$$\begin{aligned}& Z_{t}-Z_{\xi \xi }+2ZZ_{\xi }+\beta (ZY)_{\xi }=0. \end{aligned}$$
(52)

The exact solution is [47, 48]

$$\begin{aligned}& Y(\xi ,t)=a_{0}{ \bigl(1-\tanh \bigl(A(\xi -2At)\bigr) \bigr)}, \end{aligned}$$
(53)
$$\begin{aligned}& Z(\xi ,t)=a_{0}{ \biggl({ \biggl(\frac{2\beta -1}{2\alpha -1} \biggr)}- \tanh \bigl(A(\xi -2At)\bigr) \biggr)}, \end{aligned}$$
(54)

where \(a_{o}\) is an arbitrary constant, and \(A=\frac{1}{2}a_{0} (\frac{4\alpha \beta -1}{2\alpha -1} )\). The initial and boundary conditions are extracted from the exact solution. The numerical solution was computed for different values of α, β, and \(T=0.5, 1\) in the domain \([-10, 10]\) with \(a_{0}=0.05\) and \(dt=0.001\). For comparison, the error norms were computed and shown in Table 10. The stability of the scheme was discussed in terms of the spectral radius shown in the table. From the table it is clear that the proposed scheme is stable and produces better results even in large domains. The solution profile of Y and Z for \(T=1\) is plotted in Fig. 3. From the figure we can easily notice the betterment of the proposed technique.

Figure 3
figure 3

Exact and approximate solutions of Example 5 for \(T=1\), \(\alpha =0.1\), \(\beta =0.3\)

Table 10 Error norms and spectral radius for \(Y(\xi ,t)\) and \(Z(\xi ,t)\) of Example 5 for \(dt=0.001\)

Example 6

In this example, we consider one-dimensional Eq. (7) with \(\mu =-1\), \(\nu =-2\), \(\gamma =0\), and \(\alpha =\beta =1\) which leads to

$$\begin{aligned}& Y_{t}-Y_{\xi \xi }-2YY_{\xi }+(YZ)_{\xi }=0, \end{aligned}$$
(55)
$$\begin{aligned}& Z_{t}-Z_{\xi \xi }-2ZZ_{\xi }+(ZY)_{\xi }=0, \end{aligned}$$
(56)

with exact solution [22]

$$\begin{aligned}& Y(\xi ,t)=e^{-t}\sin (\xi ), \end{aligned}$$
(57)
$$\begin{aligned}& Z(\xi ,t)=e^{-t}\sin (\xi ). \end{aligned}$$
(58)

The initial and boundary conditions are extracted from the exact solution. The numerical solution was obtained for time levels \(T=0.5, 1\) in the domain \([-\pi , \pi ]\). The obtained results were compared in the form of error norms with those of the differential quadrature method [22] and are shown in Table 11. The rate of convergence using CGL collocation points is shown in Table 12. From Table 11 we can observe that the present method is stable and produces a better solution than the available techniques. In Fig. 4 the numerical and exact solutions of Y and Z are plotted for \(T=1\). The figures reflect a good agreement of the obtained numerical result with exact solution.

Figure 4
figure 4

Exact and approximate solutions of Example 6 for \(T=1\), \(\alpha =1\), \(\beta =1\)

Table 11 Error norms and spectral radius for \(Y(\xi ,t)\) and \(Z(\xi ,t)\) of Example 6 for \(dt=0.001\)
Table 12 Convergence rate of maximum error of Example 6 at \(T=1\)

Example 7

In this case, we take Eq. (7) with \(\nu =0\), \(\alpha =0\), \(\beta =0\), \(\gamma =1\), and \(\mu =-1/Re \), which leads to the following two-dimensional coupled Burger equation

$$\begin{aligned}& Y_{t}+YY_{\xi }+ZY_{\eta }-\frac{1}{Re }(Y_{\xi \xi }+Y_{\eta \eta })=0 , \end{aligned}$$
(59)
$$\begin{aligned}& Z_{t}+YZ_{\xi }+ZZ_{\eta }-\frac{1}{Re }(Z_{\xi \xi }+Z_{\eta \eta })=0 . \end{aligned}$$
(60)

The exact solution is

$$\begin{aligned}& Y(\xi , \chi , t)=0.75-0.25 \biggl[1+\exp \biggl((-4\xi +4\eta -t) \frac{Re }{32} \biggr) \biggr]^{-1}, \end{aligned}$$
(61)
$$\begin{aligned}& Z(\xi , \chi , t)=0.75+0.25 \biggl[1+\exp \biggl((-4\xi +4\eta -t) \frac{Re }{32} \biggr) \biggr]^{-1}. \end{aligned}$$
(62)

The initial and boundary conditions are taken from the exact solution. The numerical solution was computed in the domain \([0, 1]\times [0,1]\) for different values of nodal points M and Reynolds number Re when \(T=0.01\). The spectral radius and error norms were computed for \(Re =1, 10 \), and 100 and nodal points \(N=100, 400\), that is, \(N=(10 \times 10)\) and \((20 \times 20)\), and compared with the error norms obtained by Arshad [4] using the meshfree technique presented in Table 13. From the table we notice that the present results are more accurate when \(Re=1\), whereas the accuracy decreases as Re increases with increasing spectral radius. The solution and error plots are shown in Fig. 5, which shows a kink-like behavior for \(Re=50\).

Figure 5
figure 5

Solution profile and error analysis of Example 7 for \(T=0.01\), \(N=400\), \(Re =50\), and \(dt=0.005\)

Table 13 Error norms and spectral radius for \(Y(\xi , \eta , t)\) and \(Z(\xi , \eta , t)\) of Example 7 for \(T=0.01\) using regular points

Example 8

Finally, we study two-dimensional coupled Burger equations (59)–(60) in the domain \([0,1]\times [0, 1]\) with the following initial and boundary conditions taken from [4]:

$$\begin{aligned}& Y(\xi , \eta , 0)=\sin (\pi \xi )\sin (\pi \eta ), \\& Z(\xi , \eta , 0)=\bigl[\sin (\pi \xi )+\sin (2\pi \xi )\bigr] \bigl[\sin (\pi \eta )+\sin (2\pi \eta )\bigr], \\& Y(0, \eta , t)=Y(1, \eta , t)=Y(\xi , 0, t)=Y(\xi , 1, t)=0,\quad t>0, \\& Z(0, \eta , t)=Z(1, \eta , t)=Z(\xi , 0, t)=Z(\xi , 1, t)=0,\quad t>0. \end{aligned}$$

The problem was solved using the proposed technique for nodal points \(N=100, 400\), that is, \(N=(10 \times 10)\) and \((20 \times 20)\) at time \(t=0.01 \), and \(Re =1\). Due to the nonavailability of the exact solution, the obtained results were compared at different collocation points with the numerical solution by the meshfree method [4] and finite element technique [49] shown in Table 14. From the table it is clear that the proposed method produces almost the same results as those of existing methods. The solution profile for different time levels \(t=0, 0.01, 0.05\) at fixed values of η are plotted in Fig. 6. From the figure we can observe that as the time increases, \(Z(\xi , \eta , t)\)) moves from the negative part to the positive one, and the graphs tend to zero. Similarly, a 3D plot of the solution is shown in Fig. 7.

Figure 6
figure 6

Solution profile of Example 8 at \(T=0,0.01, 0.05\), \(N=400\), \(Re =1\), and \(dt=0.00125\)

Figure 7
figure 7

Solution profile of Example 8 at \(T=0.01\), \(N=400\), \(Re =1\), and \(dt=0.0001\)

Table 14 Comparison of the values of \(Y(\xi , \eta , t)\) and \(Z(\xi , \eta , t)\) of Example 8 for \(T=0.01\) using regular points

5 Conclusion

In this paper, we studied a numerical method based on the Lucas polynomials and computed solutions of three different models, including the generalized Burger–Huxley equation, generalized Burger–Fisher equation, and one- and two-dimensional nonlinear coupled Burger equations. The dependent variable is approximated by the Lucas polynomials, whereas the Fibonacci polynomials are used for its derivatives. We discussed the stability of the proposed scheme in the form of spectral radius. For comparison of the proposed method, we computed the error norms in different domains and compared the results with exact and available results in the literature. From comparison it is clear that the proposed technique gives a better accuracy and is easy to implement.