\documentclass[10pt]{article}
\usepackage{epsfig,amsfonts,amsmath,amssymb,amscd,float}
\usepackage{graphicx,caption,float,subcaption}
\usepackage{booktabs,multirow,siunitx}
\usepackage{authblk}
\usepackage{booktabs}
\usepackage{multirow}
\usepackage{soul}
\usepackage{microtype}
\usepackage{graphicx}
\usepackage{ifthen}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{xcolor}
\usepackage{epstopdf}
\DeclareMathOperator{\sech}{sech}
\pagestyle{plain} \oddsidemargin 0mm \evensidemargin 4mm \textheight 210mm
\textwidth 170mm \setlength{\baselineskip}{30pt}
\numberwithin{equation}{section}
%
\begin{document}
%
\title{Numerical study of arbitrary order Foam-Drainage equation with cubic B-spline hybrid method }
\author[1]{Alina Yousafzai\thanks{alinayousafzai1999@gmail.com}}
\author[1]{Sirajul Haq\thanks{siraj@giki.edu.pk }}
\author[2]{Abdul Ghafoor\thanks{abdulghafoor@kust.edu.pk (corresponding) }}
\author[3,4]{Kamal Shah\thanks{kamalshah@gmail.com }}
\author[4,5,6]{Thabet Abdeljawad\thanks{tabdeljawad@psu.edu.pk }}
\affil[1]{Faculty of Engineering Sciences, GIK Institute, Topi 23640, KP, Pakistan.}
\affil[2]{Institute of Numerical Sciences, Kohat University of Science and Technology,
Kohat 26000, KP, Pakistan.}
\affil[3]{Department of Mathematics and Sciences, Prince Sultan University, Riyadh, 11586, Saudi Arabia }
\affil[4]{Department of Mathematics, University of Malakand, Chakdara, Dir (L), KP, Pakistan}
\affil[5]{Department of Medical Research, China Medical Universit, Taichung40402,Taiwan}
\affil[6]{Department of Mathematics and Applied Mathematics, School of Science and Technology, Sefaco Makgatho Health Sciences University, Ga-Rankuwa, South Africa}
%\affil[6]{Department of Mathematics, Hunza Campus, Karakoram International University Gilgit, Pakistan.}
\date{}
\maketitle
\setcounter{page}{1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{abstract}
This work presents a robust and efficient numerical stratagem for the study of integer and fractional order non-linear Foam-Drainage (FD) model. The scheme first uses, usual forward difference and the $L_{1}$ formula, in integer and fractional cases, respectively. Then, the collocation approach together with cubic B-splines (CBS) basis are employed to estimate the unknown solution and its derivatives. With the help of these discretizations and Quasi-linearization, solving non-linear FD model transforms to the system of linear algebraic equations. The solution of the linear system approximates the CBS coefficients which further leads to the numerical solutions. Moreover, by Von Neumann stability it is proved that the proposed scheme is unconditionally stable. To evaluate the performance and accuracy of the technique, absolute error (AE), $L_2$, and $L_\infty$ norms are presented. The obtained outcomes are also matched with some existing results in literature. It is noted from simulations  that the proposed method gives quite accurate solutions.
\vspace{0.1cm}\\
{\bf Key words}: Non-linear Models; CBS; Finite difference approximation; Stability analysis 
\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Numerous real word problems are non-linear in nature.    
Mathematical modeling (MM) is a mechanism of representing, analyzing, and simulating these real-world systems via mathematical equations. This systematic stratagem converts the intricate physical problems into a mathematical framework, frequently in the form of equations, to predict outcomes, and understand the physical significance  of the system \cite{meerschaert2013mathematical}. Mostly, the MM leads to non-linear partial differential equations (NPDEs) of  integer and fractional order, which successfully explain a wide range of phenomena in scientific and engineering disciplines, like in control systems \cite{abbas2001using}, fluid dynamics \cite{pruess2002mathematical},  structural designs \cite{pister1972mathematical},  classical mechanics \cite{hestenes2012new}, quantum mechanics  \cite{sytnyk2021mathematical} , chemical reactions \cite{perkins2005mathematical}, molecular interactions \cite{segel1980mathematical}, material properties \cite{bhadeshia2008mathematical} in chemistry, ecosystems \cite{fulford2020mathematical} and other biological processes \cite{motta2013mathematical} etc. It also helps the environmental scientists \cite{zeidan2017mathematical} to evaluate and observe the effects of human activity on the environment.\\
In the past few years, the emergence of fractional calculus (FC) attained considerable attention to analyze the differentiation and integration to non-integer order \cite{blank1996numerical}.
In FC memory effect (ME) is the conspicuous property \cite{van2020existence},  which enhances the degree of freedom. Further, fractional derivatives are classified into global and local form, according to  ME in the system where change with respect to time occurs \cite{heydari2021numerical,podlubny1999fractional,yang2012advanced}. Recently, these derivatives have been incorporated in the modeling of various complex systems. For instance, global derivative exists in chaotic models \cite{owolabi2019formulation}, viscoelasticity \cite{yang2020general} and fluid mechanics \cite{kulish2002application}. For more details the interested readers may refer to see the related papers \cite{yu2012use,yang2017new,roohi2021chebyshev} and the reference therein. Similarly, local derivative has widespread application like, in 
 Boussinesq equation which predicts the nonlinear transverse vibration 
   mechanism \cite{ji2020fractal}, Fokker-Planck model in sub-diffusion \cite{kolwankar1998local}, Burgers' equation in audio signals propagation \cite{yang2019new}
 and fractional Kortewegde-de Vries equation  in shallow water waves \cite{yang2014modelling}.\\
Generically, the exact solutions of the NPDEs models are uncertain, and when the invoked derivatives are of arbitrary order then they become more intricate \cite{liu2022numerical}. Therefore, alternative semi-analytical and numerical approaches are essential to obtain the solutions . Some of the semi-analytical  methods include the homotopy-perturbation method (HPM) \cite{abdulaziz2008solving}, the Adomian decomposition method (ADM) \cite{jafari2006solving} and variational iteration method (VIM) \cite{momani2007numerical} etc. Similarly, various numerical strategies contain, finite difference \cite{le2016finite}, spectral \cite{rai2021spectral}, meshless \cite{belytschko1996meshless}, and  wavelet methods \cite{ghafoor2020wavelet} etc.  
Amongst the different kinfolks of numerical techniques, spline based approximations are rich  
 and is an active area of research in recent days.\\ In this work, we will analyze the following time non-linear FD model: 
     \begin{equation}\label{eq1}
     	\partial_{\tau}^{\alpha}u+2{u}^2{\partial_{x} u}- \frac{{u}{\partial_{xx}u}}{2} -{\partial_{x} u}^2=0~x\in\Gamma=[a_{ini},a_{fin}],~ \tau\geq 0,~0<\alpha\leq1,
     \end{equation}
where $\partial_{\tau}^{\alpha}u=\frac{\partial^{\alpha} u}{\partial \tau^{\alpha}},~u=u(x,\tau),~
\partial_{x} u=\frac{\partial u}{\partial x}$.  The associated boundary and initial conditions are described as follows:
\begin{eqnarray}
&& \begin{cases}
 u(a_{ini},\tau)  =\beta(\tau),\quad u(a_{fin},\tau)=\gamma(\tau), \quad \tau\geq0,  \label{eq:1.3} \\
  u(x,0)=\chi(x),~x\in[a_{ini},a_{fin}].
  \end{cases}
\end{eqnarray}       
 When  $\alpha=1$, the Eq. (\ref{eq1}) reduces to the classical FD equation. 
 In the existing studies numerous authors studied FD equation via different methods, such as ADM \cite{helal2007tanh}, homotopy analysis method \cite{singh2016homotopy}, Haar wavelet method \cite{arbabi2016semi}, reduced differential transform method (RDTM) \cite{gubes2015numerical}, meshless method \cite{shaheen2022computational} etc.	   
In numerical computations, B-spline functions particularly play a crucial role an remains and active area of research in the approximation theory and partial differential equation (PDEs). In recent past, Mittal and Arora \cite{mittal2010efficient}  developed an effective B-spline method for solving Fisher's equation . Dağ et al. \cite{daug2005numerical} provided the solution of the one-dimensional Burgers' equation using cubic B-spline (CBS). Mittal and Bhatia \cite{mittal2013numerical} proposed the B-spline method for the numerical simulations of the  second-order one-dimensional hyperbolic telegraph equation \cite{daǧ2004application}. Similarly, in \cite{kadalbajoo2012cubic} the authors explored the CBS method for the generalized Black-Scholes PDEs.\\
Inspired from the recent contribution of CBS in the numerical computations, the primal motive of this work is to put forward a new numerical strategy for the analysis and simulations of the FD equation in both integer and fractional forms. This equation has many applications in different engineering and scientific communities, like the flow of liquid  through spillways, plateau borders,  and the interrelationship of the films between the bubbles \cite{koehler2000generalized}. Moreover, the model is highly non-linear and its 
solutions is a challenging task.\\
The subsequent chunks of the paper are systematized as follows. The fundamental results related to CBS are defined in the Section \ref{sec 2}. The method description is presented in Sections \ref{sec 3} and \ref{sec 4}, respectively. In Section \ref{sec 5}, the stability of the proposed method is discussed. The test problems and conclusions are drawn in Sections \ref{sec 6} and \ref{sec 7}, respectively.

   
\section{Fundamental results}\label{sec 2}
 In this part of the manuscript, we recall the fundamentals of the CBS interpolant $\mathcal S (x)$,  which satisfy some specific end conditions. Assume that the interval  $\left[a_{ini},a_{fin}\right]$ is partitioned  into uniform sub intervals with step size $h=\frac{a_{fin}-a_{ini}}{m}$ and denoted as:
 $$\Xi=\{a_{ini}=x_{0}<x_1<x_2,\cdots,<x_m=a_{fin}\},~\mbox{where}~ x_j=a_{ini}+jh,~j=0,\cdots,m.$$
 We take the smooth CBS which is twice continuously differentiable i.e a member of the family:
   $$\tilde \Xi_{3}\Big(\Xi\Big)=\{\phi(x)|\phi(x)\in C^2[a_{ini},a_{fin}],~\mbox{where}~\phi(x)~\mbox{is a third degree poylnomial on}~ \Gamma \}.$$
 The family of CBS  $\{\mathcal B_i\left(x\right)\}_{i=-1}^{m+1}$ constitute the basis for the function defined on $\Gamma$, and are described below  \cite{splines1975variational}:
	
	\begin{equation}\label{eq5}
		{\mathcal B}_{i}(x)=
		\frac{1}{6h^3} \begin{cases}
			\left(x-x_{i-2}\right)^3 &\text{ $ \left[x_{i-2}, x_{i-1}\right]$}\\
			-3\left(x-x_{i-1}\right)^3+3h\left(x-x_{i-1}\right)^2+3h^2\left(x-x_{i-1}\right)+h^3 &\text{ $ \left[x_{i-1},x_{i}\right]$}\\
			-3\left(x_{i+1}-x\right)^3+3h\left(x_{i+1}-x\right)^2+3h^2\left(x_{i+1}-x\right)+h^3 &\text{ $ \left[x_{i},x_{i+1}\right]$}\\
			\left(x_{i+2}-x\right)^3 &\text{ $ \left[x_{i+1},x_{i+2}\right]$}\\
			0   & \text{otherwise}.
		\end{cases}
	\end{equation}
The values of $\mathcal B_{i}(x)$ and its derivatives at the nodes are shown in the following table:
	 	\begin{table}[H]
		\caption{Table of values of $\mathcal B_{i}(x)$ and its derivatives at nodes. }
\label{3T1} \vspace{-4mm}
		\begin{center}
			\begin{tabular}{c c c c c c}
				%	\multicolumn{6}{c}{$\mathrm{U}(\mathrm{w},t)$}\\
				\hline                                                            
				& $x_{i-2}$   & $x_{i-1}$ & $x_{i}$        & $x_{i+1}$        & $x_{i+2}$                           \\
				\hline
				$\mathcal B_{i}(x)$ & 0 & 1             & 4 & 1 & 0 \\[5pt]
				$\frac{d}{dx}\mathcal B_{i}(x)$ & 0 & $\frac{-3}{h}$ & 0 & $\frac{3}{h}$ & 0 \\[5pt]
				$\frac{d^2}{dx^2}\mathcal B_{i}(x)$ & 0 & $\frac{6}{h^2}$ & $\frac{-12}{h^2}$ & $\frac{6}{h^2}$ & 0 \\[5pt]
				
				\hline
			\end{tabular}
		\end{center}
	\end{table}
%Using  \eqref{eq4} and CBS \eqref{eq5}, the values of $U(x)$, $U'(x)$ and $U''(x)$ at the knots are determined in term of $c_{m}$ by 
%	\begin{align}\label{eq6}\nonumber
%		U_{m}(x) &= \frac{1}{6}(c_{m-1}+4c_{m}+c_{m+1}), \\
%		U'_{m}(x) &= \frac{1}{2h}(c_{m+1}-c_{m-1}), \\\nonumber
%		U''_{m}(x) &= \frac{1}{h^2}(c_{m-1}-2c_{m}+c_{m+1}),
%	\end{align}
%   where $U'_{m}(x)$ and $U''_{m}(x)$ denote the first and second differentiation respectively, with respect to x.
\section{Method description for integer case }\label{sec 3}
Consider Eq. (\ref {eq1}) for $\alpha=1$ as:
 \begin{equation}\label{eq7}
     	\partial_{t}u +2{u}^2{\partial_{x} u}- \frac{{u}{\partial_{xx}u}}{2} -{\partial_{x} u}^2=0,~x\in\Gamma=[a_{ini},a_{fin}],~ \tau\geq 0,~0<\alpha\leq1.
     \end{equation} 
 Using the forward approximation for time derivative, coupled with $\theta$-weighted $0\leq \theta \leq 1$ scheme,  equation \eqref{eq7} takes the following form:
 \begin{align*}\label{eq8}
 \frac{{u}_{i}^{n+1}-u^n}{d\tau}+{\theta}\left(2{u}^2{\partial_{x}u}-\frac{{u}{\partial_{xx}u}}{2}-
 {\partial_{x}u}^2\right)^{n+1}+\left(1-\theta\right)\left(2{u}^2{\partial_{x}u}-
 \frac{{u}{\partial_{xx}u}}{2}-{\partial_{x}u}^2\right)^n=0
	\end{align*}
Rearrangement of terms gives: 
	\begin{equation}\label{eq9}
		{u}^{n+1}+{\theta}{d \tau}\left(2{u}^2{\partial_{x} u}-\frac{{u}{\partial_{xx}u}}{2}-{\partial_{x}u}^2\right)^{n+1}={u}^n-\left(1-\theta\right){d \tau}\left(2{u}^2{\partial_{x}u}-\frac{{u}{\partial_{xx}u}}{2}-{\partial_{x}u}^2\right)^n
	\end{equation}
The non-linear terms $\left({u}^2{\partial_{x}u}\right)^{n+1}$, $\left({u}{\partial_{xx}u}\right)^{n+1}$ and $\left({\partial_{x}u}^2\right)^{n+1}$ are linearized using Quasi-linearization technique which are given below:
	\begin{equation}\begin{split}\label{eq10}					&\Big({u}^2{\partial_{x}u}\Big)^{n+1}=2{u}^n{\partial_{x}u}^n{u}^{n+1}
+{u}^{2n}{\partial_{x}u}^{n+1}
-2{u}^{2n}{\partial_{x}u}^n,\\ &\Big({u}{\partial_{xx}u}\Big)^{n+1}={u}^n{\partial_{x}u}^{n+1}+{u}^{n+1}{\partial_{x}u}^n
-{u}^n{\partial_{x}u}^n,\\			&\Big({\partial_{x}u}^2\Big)^{n+1}=2{\partial_{x}u}^n{\partial_{x}u}^{n+1}-{\partial_{x}u}^{2n}.\\
		\end{split}
\end{equation}
Inserting approximations from equation \eqref{eq10} in \eqref{eq9}  which gives:
	
	\begin{equation}\label{eq11}
		\begin{aligned}
			&\quad{u}^{n+1} + {\theta}{d\tau} \bigg[2\bigg(2{u}^n{\partial_{x}u}^n{u}^{n+1} + {u}^{2n}{\partial_{x}u}^{n+1}- 2{u}^{2n}{\partial_{x}u}^n\bigg) 
			-\frac{1}{2} \bigg({u}^n{\partial_{xx}u}^{n+1} + {u}^{n+1}{\partial_{xx}u}^n\\&- {u}^n{\partial_{xx}u}^n\bigg)- \bigg(2{\partial_{x}u}^n{\partial_{x}u}^{n+1} - {\partial_{x}u}^{2n}\bigg)\bigg] 
			= {u}^n - \left(1-\theta\right){d \tau} \bigg(2{u}^2{\partial_{x}u} - \frac{{u}{\partial_{xx}u}}{2} - {\partial_{x}u}^2\bigg)^n\\ 
		\end{aligned}
	\end{equation}
Further  algebraic manipulation leads to:
\begin{equation}\label{eq13}
		\begin{aligned}
			&\left(1+4\theta{d \tau}{u}^n{\partial_{x}u}^n-\frac{1}{2}\theta{d \tau}{\partial_{xx}u}^n\right){u}^{n+1}+2\theta{d \tau}\left({u}^{2n}-{\partial_{x}u}^n\right){\partial_{x}u}^{n+1}-\left(\frac{1}{2}\theta{d \tau}{u}^n\right){\partial_{xx}u}^{n+1}\\&
			={u}^n-\left(1-\theta\right){d \tau}\left(2{u}^2{\partial_{x}u}-\frac{{u}{\partial_{xx}u}}{2}-{\partial_{x}u}^2\right)^n-{\theta}{d \tau}\left(-4{u}^{2n}{\partial_{x}u}^n+\frac{{u}^n{\partial_{xx}u}^n}{2}+{\partial_{x}u}^{2n}\right).\\
		\end{aligned}
	\end{equation}
Next, using CBS,  the solution, first and second derivatives are sought in the following way:
\begin{equation}\begin{split}\label{eq14}
		&u(x,t)=\sum_{i=-1}^{m+1}{c_i(\tau)}{\mathcal {B}_i(x)},\\
        &\partial_{x}u(x,t)=\sum_{i=-1}^{m+1}{c_i(\tau)}\frac{d}{dx}{\mathcal {B}_i(x)},\\
        & \partial_{xx}u(x,t)=\sum_{i=-1}^{m+1}{c_i(\tau)}\frac{d^2}{dx^2}{\mathcal {B}_i(x)},\\
\end{split}
	\end{equation}
where $\mathcal{B}_{i}$ are the CBS and $c_i(\tau)$ stand for the time dependent coefficients to be extracted from Eq. (\ref{eq13}) with CBS collocation procedure. Substituting values from Eq. (\ref{eq14}) into Eq. (\ref{eq13}) with $x\rightarrow x_j$ one gets the following system:
\small
\begin{equation}\label{eq15}
		\begin{aligned}
&\left(1+4\theta{\Delta \tau}{u}^n{\partial_{x}u}^n-\frac{1}{2}\theta{\Delta \tau}{\partial_{xx}u}^n\right)\sum_{i=-1}^{m+1}{c^{n+1}_i(\tau)}{\mathcal{B}_i(x_j)}+2\theta{\Delta \tau}\left({u}^{2n}-{\partial_{x}u}^n\right)\sum_{i=-1}^{m+1}{c^{n+1}_i(\tau)}{\frac{d}{dx}
\mathcal{B}_i(x_j)}-\left(\frac{1}{2}\theta{\Delta\tau}{u}^n\right)\times \\&\sum_{i=-1}^{m+1}{c^{n+1}_i(\tau)}
{\frac{d^2}{dx^2}\mathcal{B}_i(x_j)}=\sum_{i=-1}^{m+1}{c^n_i(\tau)}{\mathcal{B}_i(x_j)}-\left(1-\theta\right){\Delta \tau}\left(2{u}^2\sum_{i=-1}^{m+1}{c^n_i(t)}{\frac{d}{dx}\mathcal{B}_i(x_j)}\frac{1}{2}u\sum_{i=-1}^{m+1}{c^n_i(\tau)}
{\frac{d^2}{dx^2}\mathcal{B}_i(x_j)}\right)-(1-\theta)\times \\&\Delta \tau \left(-{\sum_{i=-1}^{m+1}{c^{2n}_i(\tau)}{\frac{d}{dx}\mathcal{B}_i(x_j)}}\right)-{\theta}{\Delta \tau}\left(-4{u}^{2n}{\sum_{i=-1}^{m+1}{c^n_i(\tau)}{\frac{d}{dx}\mathcal{B}_i(x_j)}}+\frac{1}{2}{u^n}
{\sum_{i=-1}^{m+1}{c^n_i(\tau)}{\frac{d^2}{dx^2}\mathcal{B}_i(x_j)}}+{\sum_{i=-1}^{m+1}{c^{2n}_i(\tau)}
{\frac{d}{dx}\mathcal{B}_i(x_j)}}\right).\\
		\end{aligned}
	\end{equation}
The above system represents $m+1$ equations and $m+3$ unknowns. To compute the unique solution of the system, two additional conditions are required which can be adjusted using the boundary conditions as:
\begin{equation}\label{eq16}
			\left\{\begin{split}
				&\sum_{i=-1}^{m+1}c^{n+1}_i(\tau){\mathcal{B}_i(a_{ini})}=\beta(\tau),\\
				&\sum_{i=-1}^{m+1}c^{n+1}_i(\tau){\mathcal{B}_i(a_{fin})}=\alpha(\tau).
			\end{split}\right.
		\end{equation}
%\begin{equation}\begin{split}\label
%&u(x_0,\tau)=c^{n+1}_{-1}+4c^{n+1}_{0}+c^{n+1}_{1}=\beta(\tau)
%		 \Rightarrow c^{n+1}_{-1}=\beta(\tau)-\left(c^{n+1}_{1}+4c^{n+1}_{0}\right)\\& 
%u(x_m,\tau)=c^{n+1}_{N-1}+4c^{n+1}_{N}+c^{n+1}_{m+1}=\gamma(\tau)\Rightarrow
%		 c^{n+1}_{m+1}=\gamma(\tau)-\left(c^{n+1}_{N-1}+4c^{n+1}_{N}\right).
%\end{split}\end{equation}
From Eqs. (\ref{eq15}) and (\ref{eq16}) we get the following system of linear equations in compact form: 
\begin{equation}\label{eq17}
	Z_1c^{n+1}=Z_2c^{n},
\end{equation} 
where $Z_{\kappa},~\kappa=1,2$ are the square matrices with dimensions $(m+3)\times (m+3)$ and $c^{n+1},~c^{n}$ are the column matrices with dimensions $(m+3)\times 1.$ From this system, the  
unknown constants can be determined. After getting these constants, the numerical solutions can be 
estimated using Eq. (\ref{eq14}).
\section{Method description for fractional case}\label{sec 4}  
This study also focuses on time fractional problems. The fractional derivative is used in Caputo sense which is given by:
\begin{equation}\label{eq20}
		\begin{cases}
			\displaystyle	\frac{1}{\Gamma(1-\alpha)} \int_{0}^{t}
\partial_{s} u(x,s) \frac{ds}{(\tau-s)^\alpha}, &0<\alpha<1\\[7pt]
			\partial_{\tau}u, & \alpha=1.
		\end{cases}
	\end{equation} 
The Caputo derivative can be estimated via the following well-known quadrature rule \cite{zhuang2006implicit}:
\begin{equation}\label{eq:23}
		\partial_{\tau}^{\alpha}u(x,\tau ^{n+1})=\displaystyle a_\alpha\Big[{u^{n+1}-u^{n}}\Big]+a_\alpha \sum_{j=1}^{n}b_\gamma(k)\Big[{u^{n+1-j}-u^{n-j}}\Big]+\textit{O}(\delta t^{2-\alpha}),
	\end{equation}
	where
	\begin{equation*}\label{eq24}
		a_\alpha=\frac{(d \tau)^{-\alpha}}{\Gamma(2-\alpha)},~~b_\alpha(k)=\Big[(k+1)^{1-\alpha}-(k)^{1-\alpha}\Big].
	\end{equation*} 
Utilizing Eq. (\ref{eq:23}) and the strategy discussed in Section \ref{sec 3},  Eq. (\ref{eq1}) transforms to:
%Now apply this fractional time discretization to our model
%	
%	\begin{equation}\label{eq25}
%		\begin{aligned}
%			&\left(a_\alpha+4\theta{u}^n{\partial_{x}u}^n-\frac{1}{2}\theta{\partial_{xx}u}^n\right)
%\sum_{i=-1}^{m+1}{c^{n+1}_i(\tau)}{\mathcal{B}_i(x)}+2\theta\left({u}^{2n}-{\partial_{x}u}^n\right)
%\sum_{i=-1}^{m+1}{c^{n+1}_i(\tau)}{\frac{d\mathcal{B}_i(x)}{dx}}\\&-\left(\frac{1}{2}\theta{u}^n\right)
%\sum_{i=-1}^{m+1}{c^{n+1}_i(\tau)}{\frac{d^2 \mathcal{B}_i(x)}{dx^2}}		=a_\alpha\sum_{i=-1}^{m+1}{c^n_i(\tau)}{\mathcal{B}_i(x)}\\&-\left(1-\theta\right)\left(2{u}^2
%\sum_{i=-1}^{m+1}{c^n_i(\tau)}{\frac{d\mathcal{B}_i(x)}{dx}}-\frac{1}{2}u\sum_{i=-1}^{m+1}
%{c^n_i(\tau)}{\frac{d^2 \mathcal{B}_i(x)}{dx^2}}-{\sum_{i=-1}^{m+1}{c^{n}_i(\tau)}{\frac{d \mathcal{B}_i(x)}{dx}}}\right)\\&-{\theta}\left(-4{u}^{2n}{\sum_{i=-1}^{m+1}
%{c^n_i(\tau)}{\frac{d\mathcal{B}_i(x)}{dx}}}+\frac{1}{2}{u^n}{\sum_{i=-1}^{m+1}{c^n_i(\tau)}
%{\frac{d^2\mathcal{B}_i(x)}{dx^2}}}+{\sum_{i=-1}^{m+1}{c^{n}_i(\tau)}{\frac{d\mathcal{B}_i(x)}{dx}}}\right)\\&			-a_\alpha\sum_{k=1}^{n}b_\alpha(k)\Big[{u^{n+1-k}-u^{n-k}}\Big]+f^{n}_m
%		\end{aligned}
%	\end{equation}		
	\begin{equation}\label{eq26}
		\begin{aligned}
			&\sum_{i=-1}^{m+1}\Big[a_\alpha{\mathcal {B}_i(x)}+4\theta{u}^n{\partial_{x}u}^n{\mathcal{B}_i(x)}-\frac{1}{2}\theta{\partial_{xx}u}^n{\mathcal{B}_i(x)}+
2\theta\left({u}^{2n}-{\partial_{x}u}^n\right){\frac{d \mathcal{B}_i(x)}{dx}}-\frac{1}{2}\theta{u}^n{\frac{d^2 \mathcal{B}_i(x)}{dx^2}}\Big]c^{n+1}_i\\&
			=\sum_{i=-1}^{m+1}\Big[a_\alpha{\mathcal{B}_i(x)}-\left(1-\theta\right)\left(2{u}^2{\frac{d\mathcal{B}_i(x)}{dx}}-
\frac{1}{2}u{\frac{d^2 \mathcal{B}_i(x)}{dx^2}}-{\frac{d \mathcal{B}_i(x)}{dx}}^2\right)-{\theta}\left(-4{u}^{2n}{\frac{d \mathcal{B}_i(x)}{dx}}+\frac{1}{2}{u^n}{\frac{d^2 \mathcal{B}_i(x)}{dx^2}}+{\frac{d \mathcal{B}_i(x)}{dx}}^2\right)\Big]c^{n}_i\\&			-a_\alpha\sum_{k=1}^{n}b_\alpha(k)\Big[{u^{n+1-k}-u^{n-k}}\Big].
		\end{aligned}
	\end{equation}
At $x\rightarrow x_j$, the system can be written as:
	\begin{equation}\label{eq27}
		\begin{aligned} Rc^{n+1}=Tc^n-\alpha\sum_{k=1}^{n}b_\alpha(k)\Big[{u^{n+1-k}-u^{n-k}}\Big],
		\end{aligned}
	\end{equation}
where
	\begin{equation}\begin{split}\label{eq28}	&R=\sum_{i=-1}^{m+1}\Big[a_\alpha{\mathcal{B}_i(x_j)}+4\theta{u}^n{\partial_{x}u}^n{\mathcal{B}_i(x_j)}
-\frac{1}{2}\theta{\partial_{xx}u}^n{\mathcal{B}_i(x_j)}+2\theta\left({u}^{2n}-{\partial_{x}u}^n\right)
{\frac{d\mathcal{B}_{i}(x_j)}{dx}}-\frac{1}{2}\theta{u}^n{\frac{d^2 \mathcal{B}_{i}(x_j)}{dx^2}}\Big],\\			&T=\sum_{i=-1}^{m+1}\Big[a_\alpha{\mathcal{B}_i(x_j)}-\left(1-\theta\right)\left(2{u}^2{\frac{d\mathcal{B}_i(x_j)}{dx}}
-\frac{1}{2}u{\frac{d^2 \mathcal{B}_i(x_j)}{dx^2}}-{\frac{d\mathcal{B}_i(x_j)}{dx}}^2\right)-{\theta}\left(-4{u}^{2n}{\frac{d\mathcal{B}_i(x_j)}{dx}}
+\frac{1}{2}{u^n}{\frac{d^2 \mathcal{B}_i(x_j)}{dx^2}}\right)\Big]\\&-\theta{\frac{d\mathcal{B}_i(x_j)}{dx}}^2.
		\end{split}
	\end{equation}
From the above system the unknown coefficients can be extracted which can be used in Eq. (\ref{eq14}) to obtain the numerical solution.  
	\section{ Stability Analysis}\label{sec 5}
Here, stability of the proposed technique is investigated using Von Neumann stability analysis. 
Using  Eq. \eqref{eq5} and  Eq. \eqref{eq14}, the values of $u(x)$, $\partial_{x}u(x)$ and $\partial_{xx}u(x)$ at the knots in terms of $c_{i}$ are given as follows: 
	\begin{equation}\begin{split}\label{eq:a}
		u_{i}(x) &= \frac{1}{6}(c_{i-1}+4c_{i}+c_{i+1}), \\
		\partial_{x}u_{i}(x) &= \frac{1}{2h}(c_{i+1}-c_{i-1}), \\
		\partial_{xx}u_{i}(x) &= \frac{1}{h^2}(c_{i-1}-2c_{i}+c_{i+1}) .
	\end{split}
\end{equation}
For $\theta=\frac{1}{2}$ and for $\lambda=u^n,$ rewriting equation \eqref{eq7} as:
	\begin{equation}\label{eq30}
		u^{n+1}+{\Delta t}{\lambda^2}{\partial_{x}u}^{n+1}-\frac{1}{4}{\Delta \tau}\lambda{\partial_{xx}u}^{n+1}=	u^n-{\Delta \tau}{\lambda^2}{\partial_{x}u}^n+\frac{1}{4}{\Delta \tau}\lambda{\partial_{xx}u_{xx}}^n.
	\end{equation}
Substituting approximations from Eq. (\ref{eq:a}) in Eq. (\ref{eq30}) yields:
\begin{equation}\label{eq31}
		L_1c_{i-1}^{n+1}+	L_2c_{i}^{n+1}+	L_3c_{i+1}^{n+1}=L_4c_{i-1}^n+L_5c_{i}^n+L_6c_{i+1}^n,
	\end{equation}
	where
	\begin{equation}\label{eq32}
		\begin{split}
	&	L_1=\left(1-\frac{3}{h}{\Delta \tau}{\lambda^2}-\frac{3}{2h^2}{\Delta \tau}{\lambda}\right),~~~~~L_2=\left(4+\frac{3}{h^2}{\Delta \tau}{\lambda}\right),\\
	&	L_3=\left(1+\frac{3}{h}{\Delta \tau}{\lambda^2}+\frac{3}{2h^2}{\Delta \tau}{\lambda}\right),~~~~L_4=\left(1+\frac{3}{h}{\Delta \tau}{\lambda^2}+\frac{3}{2h^2}{\Delta \tau}{\lambda}\right),\\
	& L_5=\left(4-\frac{3}{h^2}{\Delta \tau}{\lambda}\right),~~~~~~~~~~~~~~~~~~~~~~~L_6=\left(1-\frac{3}{h}{\Delta \tau}{\lambda^2}-\frac{3}{2h^2}{\Delta \tau}{\lambda}\right).
		\end{split}
	\end{equation} 
Now, assume that: 
\begin{equation}\label{eq33}
		c^n=e^{ijkh}\delta^n,
	\end{equation}
where  $k$ and $h$ denote the element size and mode number, respectively, and $\delta$ is the amplification factor. From, Eqs. (\ref{eq31}-\ref{eq33}) and some algebraic manipulation one can deduce:
 \begin{equation}\label{eq34}
     		\Big[A+i B\Big]\delta^{n+1}=\Big[A-i B\Big]\delta^n,
     \end{equation}          
       where  
       	\begin{align}
       	\begin{cases}
       		\begin{split}
       			&A=2\cos(kh)+\left(4+\frac{3}{h^2}{\Delta t}{\lambda}\right),\\
       			&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\
       			&B=\left(\frac{6}{h}{\Delta t}{\lambda^2}+\frac{3}{h^2}{\Delta t}{\lambda}\right)\sin(kh),
       		\end{split}
       	\end{cases}
       \end{align}
       and
       \begin{equation}
       \left| \frac{\delta^{n+1}}{\delta^n}\right| =\frac{|A-i B|}{|A+i B|}.
       \end{equation}
From the above statement we can write
        \begin{equation}
       	|\delta| =\sqrt{\frac{A^2+B^2}{A^2+B^2}}=1,
       \end{equation}
which shows that the scheme is unconditionally stable. 
\section{Illustrative examples}\label{sec 6}
Here, the proposed scheme is implemented for the solutions of some test problems. 
To measure the accuracy of the technique, the following error norms are employed:  
	\begin{itemize}
		\item Infinity error norm: 
		\begin{equation*}\label{eq35}
			L_{\infty}=\max_{1\leq i\leq m} |u_E(\mathrm{x_i},\mathrm{t_n})-u_A(\mathrm{x_i},\mathrm{t_n})|,
		\end{equation*}
		\item $L_2$ error norm:
		\begin{equation*}\label{eq36}
			L_2=\sqrt{h \sum_{i=1}^{m}(u_E(\mathrm{x_i},\mathrm{t_n})-u_A(\mathrm{x_i},\mathrm{t_n}))^2},
		\end{equation*}
		\item Absolute error:
		\begin{equation*}\label{eq37}
			AE=|u_E(\mathrm{x_i},\mathrm{t_n})-u_A(\mathrm{x_i},\mathrm{t_n})|,
		\end{equation*}
	\end{itemize}
where $u_E$ and $u_A$ denote the exact and approximate solutions, respectively.
\subsection*{Example 1}\label{E1}
Consider Eq. \eqref{eq1} in integer form with exact solution as follows \cite{behzadi2014use}:
\begin{equation}\label{eq38}
	\begin{split}
	&u(x,t)=-\sqrt{k}\tanh\left(\sqrt{k}\left(x-k\tau\right)\right),
	\end{split}
\end{equation}
where $k$ represents the speed of the wave front.
The corresponding initial and boundary conditions are extracted from the given exact solution.
The problem is solved via the proposed technique and the results are recorded. In Table \ref{T2}
the absolute error at various points are reported for $k=1$ and $t=0.001$. Similarly, in Table \ref{T3} the absolute error is noted for $k=2$ and $\tau=0.0001$. From first table it is visible that the results of the proposed method are comparable with the available results obtained through
radial basis functions method and hybrid method \cite{shaheen2022computational,arslan2020numerical}, and in second table our results are better.
In Fig. \ref{1f1} the three dimensional solutions profiles with absolute error are displayed for $k=,~\tau=0.001.$ Also, in Fig.  \ref{3f1} the two  dimensional solution profiles are presented. From figures it is clear that computed solution mutually agrees with the exact solution. 
		\begin{table}[H]
		\caption{Results comparison of example 1 at $\tau=0.001$ and $k=1$. }\label{2T1} \vspace{-4mm}
		\begin{center}
			\begin{tabular}{c c c c c c}\hline                                                            
				$x$     & Present method   & \multicolumn{1}{c}{HM\cite{arslan2020numerical}} & IQ\cite{shaheen2022computational}       & \multicolumn{1}{c}{MQ \cite{shaheen2022computational}}        & GA \cite{shaheen2022computational}                            \\
				\hline
			0 &	1.00E-04 &	5.79950E-04&	8.98874E-04&	9.00000E-04&	6.38429E-04 \\
			0.2 &	4.51E-04 &	8.20360E-04&	7.62965E-04&	7.33814E-04&	5.33519E-04 \\
			0.4 &3.36E-04 &	9.38980E-04&	5.25293E-04&	4.09539E-04&	3.77066E-04 \\
			0.6 &1.46E-04 &	9.09350E-04	&2.93868E-04&	7.89034E-05&	6.73231E-06 \\
			0.8 &1.09E-05&	7.82290E-04	&2.89867E-06&	1.78086E-05&	2.09170E-04 \\
			1 & 4.21E-05&	6.40590E-04&	3.52396E-04&	3.78294E-04&	4.51210E-05 \\
			
				
				\hline
			\end{tabular}
\label{T2}
		\end{center}
	\end{table}
	
		\begin{table}[H]
		\caption{Results comparison of example 1 at $\tau=0.0001 $ and $k=2$. }\label{3T1} \vspace{-4mm}
		\begin{center}
			\begin{tabular}{c c c c c c}\hline                                                            
				$x$     & Present method   & \multicolumn{1}{c}{HM\cite{arslan2020numerical}} & IQ\cite{shaheen2022computational}       & \multicolumn{1}{c}{MQ \cite{shaheen2022computational}}        & GA \cite{shaheen2022computational}                            \\
				\hline
			0&	4.00E-05&	1.39360E-04&	3.49451E-04&	3.60000E-04&	3.58354E-04 \\
			0.2&	9.58E-05&	3.41610E-04&	1.60253E-04&	2.50909E-04&	2.65987E-04 \\
			0.4&	1.29E-05&	4.37679E-04&	1.44167E-04&	6.29466E-05&	1.18012E-04 \\
			0.6&	4.95E-05&	4.20165E-04&	9.59463E-05&	6.39323E-05&	2.99209E-06 \\
			0.8&	6.66E-05&	3.45387E-04&	8.30372E-05&	1.01200E-04&	2.91812E-05 \\
			1&	8.43E-06&	2.61457E-04&	5.90165E-05&	7.58986E-05&	7.51460E-05 \\
				
				\hline
			\end{tabular}
\label{T3}
		\end{center}
	\end{table}
	
		\begin{figure}[H]
		\begin{center}
			\begin{subfigure}[H]{0.45\textwidth}
				\centering
				\includegraphics[width=2.6in, height=1.5in]{SA1.eps}\vspace{0mm}
				\caption{Approximate}\label{af1}
			\end{subfigure}
			\begin{subfigure}[H]{0.45\textwidth}
				\centering
				\includegraphics[width= 2.6in, height=1.5in]{SE1.eps}\vspace{0mm}
				\caption{Exact}\label{bf1}
			\end{subfigure}
				\begin{subfigure}[H]{0.45\textwidth}
				\centering
				\includegraphics[width=2.6in, height=1.5in]{E1.eps}\vspace{0mm}
				\caption{AE}\label{cf1}
			\end{subfigure}
				\caption{Three dimensional solution profiles of example 1 with absolute error. }\label{1f1}	
		\end{center}	
	\end{figure}
	
	\begin{figure}[H] 
		\begin{center}
				\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width=2.6in, height=1.5in]{EVA11.eps}\vspace{0mm}
					\end{subfigure}
				\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width=2.6in, height=1.5in]{EVA12.eps}\vspace{0mm}
					\end{subfigure}
				\caption{Exact versus approximate solution for $k=2$ and $\tau=0.0001.$}\label{3f1}
			\end{center}
	\end{figure}
	\subsection*{Example 2}
	Consider Eq. \eqref{eq1}, in integer form with the following exact solution:
	\begin{equation}\label{eq39}
		\begin{split}
			&u(x,t)=\frac{-1}{2}+\frac{1}{1+e^{x-\frac{\tau}{4}}}.
		\end{split}
	\end{equation}
 Simulations are computed for different parameters. In Table \ref{T4} the attained
 absolute error at distinct points are presented. Likewise, the computed error norms for 
 different values of $m=2^j,~j=1,\cdots,5,$ and various times are shown in Table \ref{T5}.
 Simulations indicate that proposed method works well for the given problem. Moreover, the 
 exact and approximate measures are plotted in Fig. \ref{3f3} together with an absolute error, 
 in both two and three dimensional forms, respectively. Graphical view discloses that exact and approximate solutions are promises well.      
\begin{table}[H]
		\caption{Absolute error for $h=0.1$ and $\tau=0.001$ of example 2.}\label{3T4}
		\begin{center}
			\begin{tabular}{c c c c c c}\hline
			
				$x$  & Exact  & Approximate  &  AE    \\\hline
				0  & 0.00006250 &0.00005625  & 6.2500E-06 \\
				0.2  & -0.0498 & -0.0499 & 1.3701E-04 \\
				0.4& -0.0986 & -0.0987 & 1.2013E-04 \\
				0.6& -0.1456 & -0.1457 & 1.0961E-04 \\
				0.8& -0.1899 & -0.1900 & 9.6098E-05\\
				1& -0.2310 & -0.2310 & 4.9158E-06\\
				\hline
			\end{tabular}\label{T4}
		\end{center}
	\end{table} 
\begin{table}[H]
		\caption{Error norms of example 2 for different parameters and $\tau=0.01$.}\label{3T5}
		\begin{center}
			\begin{tabular}{c c c c c c} \hline \hline
			&	\multicolumn{2}{c}{$\Delta \tau=0.001$ } & \multicolumn{2}{c}{ $\Delta \tau=0.0001$}                 \\
				\hline
				$m$  & $L_2$  & $L_\infty$  & $L_2$  & $L_\infty$     \\\hline
				2  & 7.09E-04 & 2.00E-02 & 6.86E-04 &9.70E-04 \\
				4  & 9.67E-04 & 1.20E-03 & 9.62E-04 &1.30E-03\\
				6  & 1.00E-03 &  1.30E-03 & 1.00E-03 &1.40E-03\\
				8 & 1.10E-03 & 1.40E-03 & 1.10E-03 &1.40E-03\\
				10 & 1.10E-03 & 1.40E-03 & 1.10E-03&1.40E-03\\
				\hline
			\end{tabular}\label{T5}
		\end{center}
	\end{table}
\begin{figure}[H] 
		\begin{center}
				\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width=2.6in, height=1.5in]{SE2.eps}\vspace{0mm}
						\caption{Exact}
					\end{subfigure}
				\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width= 2.6in, height=1.5in]{SA2.eps}\vspace{0mm}
						\caption{Approximate}
					\end{subfigure}
				\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width= 2.6in, height=1.5in]{EVA2.eps}\vspace{0mm}
						\caption{Exact versus approximate}
					\end{subfigure}
						\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width= 2.6in, height=1.5in]{E2.eps}\vspace{0mm}
						\caption{AE}
					\end{subfigure}
				\caption{Exact versus approximate solutions of example 2 absolute error.}\label{3f3}
			\end{center}
	\end{figure}
\subsection*{Example 3}
In this example, fractional form of Eq. \eqref{eq1} is considered. The exact solution for fractional FD equation has been taken from \cite{habib2021numerical} which is given below:
\begin{equation}\label{eq40}
	u(x,t)=-\sqrt{k}\tanh({\sqrt{k}x})+k^2\sech^2(\sqrt{k}x)
\left[\frac{\tau^\alpha}{\Gamma(1+\alpha)}\right]+k^{\frac{7}{2}}\sech^2(\sqrt{k}x)\tanh
(\sqrt{k}x)\left[\frac{\tau^\alpha}{\Gamma(1+\alpha)}\right]^2.
	\end{equation}
All necessary conditions are used from the analytical solution. Simulation are performed for various values of fractional order derivative. In Tables \ref{T6}-\ref{T8}, the compute solutions along with absolute error for $(\alpha=0.7,0.8,0.9)$ at different mesh points are matched with the work presented in \cite{shaheen2022computational}. It is noted that proposed method performs well. 
The efficiency of the method is also portrayed using the graphical representation of the solutions. In Fig. \ref{F4} the two dimensional solutions pattern and its corresponding error are presented. The three dimensional solutions pattern for $\alpha=0.8,0.9$ are displayed in Fig. \ref{f5}. From figures it is obvious that calculated solution are pretty much closed to exact solution. The simulations of this problem also show that the suggested scheme is quite suitable and easily implementable for such fractional models.               
\begin{table}[H]
		\caption{Results comparison of example 3 $\alpha=0.7$, $\tau=0.01.$ }\label{3T6}
		\begin{center}
			\hspace*{-0.72cm}	\begin{tabular}{c c c c c c} \hline \hline
			%	\multicolumn{4}{c}{ }                           \\
				$x$& Exact & Approximate & Present method& IQ \cite{shaheen2022computational} & MQ \cite{shaheen2022computational} \\ 	\hline
				-4&	1.732048652&	1.732048591&	6.04E-08&	1.80E-04&	1.39E-02 \\
				-2&	1.729852038	&1.730436162	&5.84E-04&	6.59E-04&	7.55E-04 \\
				0&	0.394322232&	0.038939663	&3.55E-01&	7.41E-02&	0.407953 \\
				2&	-1.726767357&	-1.732521165&5.75E-03&	9.53E-03&	8.93E-03 \\
				4&	-1.732045623&	-1.732045778&	1.55E-07&	1.54E-03&	1.65E-04 \\
				
				\hline
			\end{tabular}\label{T6}
		\end{center}
	\end{table}
	\begin{table}[H]
		\caption{Results comparison of example 3 at $\alpha=0.8$, $\tau=0.01.$}\label{3T7}
		\begin{center}
			\begin{tabular}{c c c c c c c c } \hline \hline
			%	\multicolumn{4}{c}{$t=0.01, \alpha=0.8$}                           \\
		
			$x$& Exact & Approximate & Present method& IQ \cite{shaheen2022computational} & MQ \cite{shaheen2022computational} \\ 	\hline
			-4&	1.73204828&	1.73204822&	5.51E-08&	8.23E-05&	1.71E-04 \\
			-2&	1.72947674&	1.72997423&	4.97E-04&	1.75E-04&	3.95E-04 \\
			0&	0.24272462&	0.20755948&	3.51E-02&	1.21E-01&	1.33E-01 \\
			2&	-1.72757797&	-1.72232369&	5.25E-03&	1.12E-03&	7.58E-04 \\
			4&	-1.73204642&	-1.73204651&	9.56E-08&	4.31E-04&	7.01E-05 \\
			\hline	
			\end{tabular}\label{T7}
		\end{center}
	\end{table}
	
		\begin{table}[H]
		\caption{Results comparison of example 3 $\alpha=0.9$, $\tau=0.01.$ }\label{3T8}
		\begin{center}
			\begin{tabular}{c c c c c c c c } \hline \hline
			%	\multicolumn{4}{c}{$t=0.01, \alpha=0.9$}                           \\
				
				$x$& Exact & Approximate & Present method& IQ \cite{shaheen2022computational} & MQ \cite{shaheen2022computational} \\ \hline
				-4&	1.73204755&	1.73204754&	6.34E-09&	2.78E-06&	7.08E-08 \\
				-2&	1.72873238&	1.72879112&	5.87E-05&	1.15E-04&	1.10E-04 \\
				0&	0.01867124&	0.01529607&    	3.37E-03&	1.18E-01&	1.15E-01 \\
				2&	-1.72858632&	-1.72852372&	6.26E-05&	2.57E-05&	6.26E-07 \\
				4&	-1.73204741&	-1.73204741&	6.62E-09&	4.24E-05&	4.82E-05 \\
				
				\hline	
			\end{tabular}\label{T8}
		\end{center}
	\end{table}
\begin{figure}[H] 
		\begin{center}
			\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width=2.6in, height=1.5in]{EVA3.eps}\vspace{0mm}
						\caption{Exact and approximate solutions at $\alpha=0.9$}		
					\end{subfigure}	
						\begin{subfigure}[H]{0.45\textwidth}
						\centering
						\includegraphics[width=2.6in, height=1.5in]{{AE32.eps}}\vspace{0mm}
						\caption{AE at $\alpha=0.9$}	
					\end{subfigure}	
					\caption{Exact and approximate solutions with corresponding error of example 3.} \label{F4}
			\end{center}
	\end{figure}

\begin{figure}[H] 
		\begin{center}
			\begin{subfigure}[H]{0.45\textwidth}
				\centering
				\includegraphics[width=2.6in, height=1.5in]{SE331.eps}\vspace{0mm}
				\caption{Exact solution for $\alpha=0.8$}
			\end{subfigure}
			\begin{subfigure}[H]{0.45\textwidth}
				\centering
				\includegraphics[width= 2.6in, height=1.5in]{SA331.eps}\vspace{0mm}
				\caption{Approximate solution for $\alpha=0.8$}
			\end{subfigure}
			\begin{subfigure}[H]{0.45\textwidth}
				\centering
				\includegraphics[width= 2.6in, height=1.5in]{SE332.eps}\vspace{0mm}
				\caption{Exact solution for $\alpha=0.9$}
			\end{subfigure}
			\begin{subfigure}[H]{0.45\textwidth}
				\centering
				\includegraphics[width= 2.6in, height=1.5in]{SA332.eps}\vspace{0mm}
				\caption{Approximate solution for $\alpha=0.9$}
			\end{subfigure}
			\caption{Three dimensional graphical solutions of example 3.}\label{f5}
		\end{center}
	\end{figure}
\section{Conclusions}\label{sec 7} 
In this paper a hybrid numerical method, consisting of CBS and finite difference approximations, has been proposed for the solution of FD  equation. The key features of the technique includes, the formation of the linear system of equations, unconditional stability and quite accurate solutions. Three test problems have been solved with the help of the proposed method. The effectiveness of the scheme has been judged via using two error measures namely: $L_{2}$ norm and $L_{\infty}$ norm  and absolute error.
 The obtained results have been compared with other numerical techniques in existing studies. Further, the simulation results have been presented graphically. From tabulated recorded data and different solution  profiles it is stated the proposed works well for the aforementioned models and would be useful to solve such other systems in different engineering systems.      
\section*{Declarations}
\subsection*{Competing interests}
The authors declare that they have no conflict of interest.
\subsection*{Funding}
\section*{Funding}
This project is funded by Prince Sultan University, Riyadh, Saudi Arabia.
\subsection*{Data Availability}
The data would be provided on demand.
\newpage
\bibliography{Bibtex}
\bibliographystyle{ieeetr}
\end{document}
%%%%%%%%%%% Table-3 problem 1 %%%%%%%%%%%%%%
