Spectral Optimized Multiderivative Hybrid Block Method for Fitzhugh–Nagumo Equations

Abstract

The Fitzhugh–Nagumo equation, a key model for excitable systems in biology and neuroscience, requires efficient numerical methods due to its nonlinear nature. A spectral optimized multiderivative hybrid block method is proposed, constructed using a multistep collocation and interpolation technique with an approximated power series as the basis function. Incorporating two optimal intra-step points, the method demonstrates improved accuracy, with its consistency, convergence, and absolute stability rigorously analyzed. By combining the optimized multiderivative hybrid block method in time with a spectral collocation method in space, the approach demonstrates potency and flexibility in solving partial differential equations. Prior to using the spectral method, the partial differential equation is linearized using a linear partition technique. Numerical experiments confirm the accuracy and efficiency of the method compared to existing methods, demonstrating the potential of the method as a robust framework for solving partial differential equations requiring both high accuracy and stability.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Nonlinear PDEs are extensively utilized for describing complex physical occurrences in diverse scientific fields, including fluid mechanics. These nonlinear models have also found importance in areas like biology and economics (see [1, 2]). One well-known example is the Fitzhugh–Nagumo (F-N) equation, which was developed by previous researchers [3, 4].$$\begin{array}{}\text{(1)}& w_t=w_{xx}+w1-ww-\delta ,\end{array}$$where $0\le \delta \le 1$ and $wx,t$ denotes the dependent function that relies on both the time variable $t$ and the space variable $x$. Equation (1) has various applications in physics, applied mathematics and mathematical modeling of complex physical and biological systems (see [5–8]). This equation incorporates both the diffusion and non-linearity, which is governed by the term $w1-ww-\delta$.

The F-N equation with diffusion arose as a special case of the famed Hodgkin–Huxley model, which earned the Nobel Prize in 1963. Nagumo’s work [4] paved the way for extensive research in neuroscience, leading to the development of various extended forms of the F-N model. This system is widely used to describe how a voltage of neuron changes over time and space and is given by$$\begin{array}{c}\text{(2)}& w_t=\mathrm{\Lambda }{w}_{xx}-y+w-\frac{w^3}{3},y_t=\epsilon -\mathrm{\Delta }y+w,\end{array}$$where $\mathrm{\Lambda },\epsilon ,$ and $\mathrm{\Delta }$ are real parameters. The functions $wx,t$ and $yx,t$, represent the voltage potential (fast variable) and the recovery variable (slow variable), respectively. This model is well-established in neurodynamics [9, 10]. An exact solution to the F-N model (2) was presented in [9].

Several methods have been investigated for the direct solution of equation (1), such as the adomian decomposition method (ADM) [11], variational iteration method (VIM) [12, 13], shifted Legendre polynomial (SLP) [14] and differential transform method (DTM) [15], including the Haar wavelet method presented in [16]. According to [17], analytical solutions for PDEs related to Nagumo reaction-diffusion have been derived using the variational method. Additionally, the homotopy analysis technique has effectively been utilized to provide approximate solutions for the standard F-N equation [18]. Furthermore, the approximate conditional symmetry method as proposed in [19], has been utilized to find approximate solutions for the perturbation of equation (1). Exact solutions have also been obtained using the first integral method [20] and Jacobi elliptic function [21]. While some of these methods require many terms for accurate analytical solutions, others face challenges in determining a suitable initial guess.

Among the diverse array of approaches used to solve PDEs (see [22–24]), the nonstandard finite difference (NSFD) techniques have recently emerged as one of the most effective approaches [25]. The exact finite-difference method (see [26–29]) is a special NSFD scheme. These methods produce reliable results across a wide range of initial and boundary conditions (BCs). However, they can be computationally demanding and prone to numerical instability.

Spectral methods are a class of numerical techniques that approximate solutions to differential equations (DEs). The three traditional spectral methods are the Galerkin, Tau and collocation methods [30]. Spectral techniques have attracted interest due to their computational efficiency and the relative simplicity with which they are used to solve nonlinear DEs [31, 32]. These methods are based on expanding the solution to the equation in terms of a set of basis functions like polynomial kind or Fourier basis [33]. Spectral collocation methods (SCMs), such as the Lagrange pseudo-spectral approach [34], the uniform Haar wavelets and quasilinearization process [35], and the Chebyshev SCMs (see [36–38]), have been presented for solving time-dependent DEs, including nonlinear and stochastic Burgers’ equation. The reaction-diffusion Brusselator system and several nonlinear Burgers’ type equations were effectively solved numerically utilizing the differential quadrature method (DQM) [39, 40]. Motsa et al. [41] introduced the bivariate SCM with the use of Chebyshev and Lagrange cardinal polynomials. This approach involves applying SCM to the time variable to transform the PDE into an ODE, followed by the spatial direction to form systems of linear equations and simplifying the solution process. Initially applied to evolutionary parabolic PDEs, this method has been adapted in [42] to solve various classes of PDEs with two independent variables. The precision of SCMs relies on the choice of grid points and the type of basis functions used to construct the interpolating polynomial [43]. However, their effectiveness decreases when applied to PDEs over extended time intervals. By using fewer grid points, SCMs not only enhance accuracy but also significantly reduce computational costs [44]. The ongoing pursuit of greater precision and computational efficiency continues to drive the development of new techniques. Spectral methods have been combined with other numerical techniques to enhance performance when utilized to solve various types of DEs. One such approach, presented in [45], integrates the finite difference method for time discretization with the spectral method for spatial discretization to solve DEs. Despite the advancements, existing methods often struggle to balance high accuracy, stability, and computational efficiency for nonlinear PDEs like the F-N equation.

Hybrid block method (HBM) has been utilized to solve equation (1) by transforming it into an equivalent system of ordinary DEs through the method of lines [46]. It was noted that as the value of time increases, the error also increases, indicating that shorter time is required to achieve more accurate solution approximations. The Hermite block method introduced in [47] was developed using Hermite polynomials as basis functions to solve second-order nonlinear elliptic PDEs, aiming to enhance numerical accuracy and efficiency. Olaoluwa Omole et al. [48] employed an implicit block method to solve PDEs, achieving greater precision and faster convergence.

Another approach is the multiderivative (MDHBMs), which offer a significant advancement in computational techniques for solving DEs. These methods incorporate higher-order derivatives in the formulation to enhance accuracy and stability without significantly increasing computational complexity [49–51]. By utilizing second or higher-order derivative terms, multiderivative methods reduce truncation errors and improve the convergence rates [51–53]. By merging the strengths of linear multistep methods and block techniques, these methods can efficiently solve both linear and nonlinear problems. Kayode $\text{et\hspace{0.17em}al}.$ [52] introduced a predictor-corrector technique utilizing higher derivative to solve IVPs. To achieve optimal accuracy in block hybrid methods (BHMs), various optimized BHMs for the numerical integration of IVPs have been presented in [54–57]. A one-step optimized HBM for solving systems of IVPs was presented in [58]. An efficient optimized HBM for integrating second-order IVPs was presented in [59]. Recently, Yakubu and Sibanda [60] used a power series polynomial as a basis function to develop an optimized second-derivative HBM to maximize accuracy when utilized to solve IVPs. While these methods have shown promise for initial value problems (IVPs), their application to nonlinear PDEs have not been extensively studied. This gap highlights the need for innovative numerical techniques to tackle nonlinear PDEs with greater efficiency and accuracy.

To address these challenges, this study introduces a spectral optimized MDHBM (SOMDHBM) for solving the F-N equations (1) and (2). By combining the high accuracy of spectral methods with the stability and computational efficiency of multiderivative hybrid block techniques, the SOMDHBM offers a robust framework for solving nonlinear and diffusive characteristics of the equation. The F-N equation is first linearized using a linear partition mapping technique before applying the SOMDHBM. The effectiveness of the SOMDHBM is validated through numerical experiments, demonstrating the accuracy and efficiency compared to existing techniques such as the polynomial DQM (PDQM) [61], compact finite difference block method (CFDBM) [49], ADM [11], VIM [12, 13], DTM [15], exact finite-difference and NSFD scheme [62], İnan $\text{et\hspace{0.17em}al}.$ [63] and exponential finite difference method (ExpFDM) [64].

2. Derivation of the Method

In this section, the SOMDHBM is developed for a parabolic PDE of the form$$\begin{array}{}\text{(3)}& \dot{w}=fx,t,w,w^{\prime },w^{″},\end{array}$$with the initial and BCs$$\begin{array}{}\text{(4)}& wx,0=w_{0}x,wa,t={\aleph }_{a}t,wb,t={\aleph }_{b}t,\end{array}$$where ${\aleph }_{a}t$ and ${\aleph }_{b}t$ are time-dependent functions or constants. The dot and the prime denote partial differentiation with respect to $t$ and $x$, respectively. The SOMDHBM is derived from the combination of OMDHBM and the SCM, with an appropriate linear transformation to convert $x\in a,b$ to $\widehat{x}\in -1,1$. We develop the OMDHBM for the DEs of the form$$\begin{array}{}\text{(5)}& \dot{w}t=ft,wt,w{t}_0=w_0,\end{array}$$where the dot denotes the first-derivative of $w$ with respect to the time $t$. The solution $wt$ of equation (5) is approximated at four collocation points$$\begin{array}{}\text{(6)}& t_{n+p_i}=t_n+h{p}_i,\end{array}$$in the closed interval $t_n,t_{n+1}$, $h$ is the chosen step-size, $p_i$ is the collocation parameters, $i=0,1,\dots ,3$ and $n$ is the grid index. The following points will be utilized$$\begin{array}{}\text{(7)}& t_{n+p_0},t_{n+p_1},t_{n+p_2},t_{n+p_3},\end{array}$$where $p_0=0$ and $p_3=1$. The points $p_1=c$ and $p_2=1-c$ are inter-step points. These choices will ensure that the method is stable and achieves higher accuracy. The solution $wt$ will be approximated by a power series polynomial of the form$$\begin{array}{}\text{(8)}& wt\approx Wt={\displaystyle \sum }_{j=0}^5{\eta }_{n,j}{t-t_n}^j,\end{array}$$where the unknown coefficients ${\eta }_{n,j}$ is obtained through a system of five equations with five unknowns generated from the equations$$\begin{array}{c}\text{(9)}& {\dot{W}}_{n+p_i}={\displaystyle \sum }_{j=1}^{5}j{\eta }_{n,j}{t-t_n}^{j-1}=f{t}_{n+p_i},w_{n+p_i},i=0,1,\dots ,3,{\ddot{W}}_{n+p_3}={\displaystyle \sum }_{j=2}^{5}jj-1{\eta }_{n,j}{t-t_n}^{j-2}=g{t}_{n+p_3},{\dot{w}}_{n+p_3},\end{array}$$with the initial condition$$\begin{array}{}\text{(10)}& W{t}_0={\eta }_{0,j}=w_0.\end{array}$$

The Mathematica 13.0 code presented in Algorithm 1, for equations (8)–(10), is used to determine the unknown coefficients ${\eta }_{n,j}$.

Algorithm 1: Mathematica code for obtaining the coefficients ${\mathit{\eta }}_{\mathit{n},\mathit{j}}$.

1. Define Points: ${\mathit{p}}_{\mathit{i}}=\mathbf{0},\mathit{c},\mathbf{1}-\mathit{c},\mathbf{1}$; Set: $\mathit{J}=\mathbf{\text{Length}}\mathit{p}$;

2. Collocate: ${\mathit{t}}_{{\mathit{p}}_{\mathit{i}}+\mathit{n}}=\mathit{h}{\mathit{p}}_{\mathit{i}}+\mathit{n}$, for $\mathit{i}=\mathbf{0},\dots ,\mathit{J}$;

3. Set:

  $\mathbf{W}=\mathbf{\text{Sum}}{\mathit{\eta }}_{\mathit{n},\mathit{j}}{\mathit{t}-{\mathit{t}}_{\mathit{n}}}^{\mathit{j}},\mathit{j},\mathbf{0},\mathit{J}+\mathbf{1}$;

  ${\mathbf{W}}^{\prime }=\mathbf{\text{Table}}\mathbf{\text{Sum}}\mathit{j}\mathbf{·}{\mathit{\eta }}_{\mathit{n},\mathit{j}}\mathbf{·}{\mathit{t}-{\mathit{t}}_{\mathit{n}}}^{\mathit{j}-\mathbf{1}}/.\mathit{t}⟶{\mathit{t}}_{\mathit{n}+{\mathit{p}}_{\mathit{i}}},\mathit{j},\mathbf{1},\mathit{J}+\mathbf{1}=={\mathit{f}}_{\mathit{n}+{\mathit{p}}_{\mathit{i}}}$, for $\mathit{i}=\mathbf{1},\dots ,\mathit{J}$;

  ${\mathbf{W}}^{″}=\mathbf{\text{Table}}\mathbf{\text{Sum}}\mathit{j}\mathbf{·}\mathit{j}-\mathbf{1}\mathbf{·}{\mathit{\eta }}_{\mathit{n},\mathit{j}}\mathbf{·}{\mathit{t}-{\mathit{t}}_{\mathit{n}}}^{\mathit{j}-\mathbf{2}}/.\mathit{t}⟶{\mathit{t}}_{\mathit{n}+{\mathit{p}}_{\mathit{i}}},\mathit{j},\mathbf{2},\mathit{J}+\mathbf{1}=={\mathit{g}}_{\mathit{n}+{\mathit{p}}_{\mathit{i}}}$, for $\mathit{i}=\mathit{J}$;

4. Identify Unknowns: $\mathbf{\text{Unknowns}}=\mathbf{\text{Table}}{\mathit{\eta }}_{\mathit{n},\mathit{i}},\mathit{i},\mathbf{0},\mathit{J}+\mathbf{1}$;

5. Apply Initial Condition: $\mathbf{\text{Initial}}=\mathit{W}/.\mathit{t}⟶{\mathit{t}}_{\mathit{n}}=={\mathit{w}}_{\mathit{n}}$;

6. Combine All Equations: $\mathbf{\text{AllEquations}}=\mathbf{\text{Join}}\mathbf{W}’,\mathbf{W}”,\mathbf{\text{Initial}}$;

7. Solve for Unknowns:

  $\mathbf{\text{Solutions}}=\mathbf{\text{Simplify}}\mathbf{\text{Flatten}}\mathbf{\text{Solve}}\mathbf{\text{AllEquations}},\mathbf{\text{Unknowns}}$;

Substituting ${\eta }_{n,j}$ in the continuous approximation (8) and evaluating at the collocation points (6). The method obtained is given below:$$\begin{array}{c}\text{(11)}& w_{n+c}=\frac{hc}{60{c-1}^{2}2c-1}-3c^3{f}_{n-c+1}+27c^3{f}_{n+c}+13c^2{f}_{n-c+1}-85c^2{f}_{n+c}+16c^5-94c^4+233c^3-275c^2+150c-30f_n\\ +-16c^5+94c^4-137c^3+47c^2+20c-10f_{n+1}-20c{f}_{n-c+1}+90c{f}_{n+c}\\ +10f_{n-c+1}-30f_{n+c}+\frac{h^2{c}^{2}4c^2-10r+5g_{n+1}}{30c-1}+w_n,\\ \text{(12)}& w_{n+1-c}=-\frac{hc-1}{60c^{2}2c-1}27c^3{f}_{n-c+1}-3c^3{f}_{n+c}+4c^2{f}_{n-c+1}-4c^2{f}_{n+c}+c16c^4+14c^3+17c^2-20c+3f_n\\ -16c^5+14c^4-79c^3+40c^2+c-2f_{n+1}+c{f}_{n-c+1}-3c{f}_{n+c}-2f_{n-c+1}\\ -\frac{h^{2}4c^2+2c-1{c-1}^2{g}_{n+1}}{30c}+w_n,\\ \text{(13)}& w_{n+1}=\frac{h}{60{c-1}^2{c}^{2}2c-1}-5c^2{f}_{n-c+1}+5c^2{f}_{n+c}+40c^4-100c^3+86c^2-29c+3c{f}_n\\ +80c^5-200c^4+154c^3-31c^2-7c+2f_{n+1}+7c{f}_{n-c+1}-3c{f}_{n+c}-2f_{n-c+1}\\ -\frac{h^{2}5c^2-5c+1g_{n+1}}{30c-1c}+w_n.\end{array}$$

The method can also be presented in the form$$\begin{array}{}\text{(14)}& w_{n+p_i}=w_n+h{\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{f}_{n+p_i}+h{\displaystyle \sum }_{j=3}^1{\beta }_{i,j}{g}_{n+p_j},i=1,2,3,\end{array}$$and in matrix form$$\begin{array}{}\text{(15)}& E_1{W}_{n+1}=E_0{W}_n+h{A}_1{F}_{n+1}+A_0{F}_n+h{B}_1{G}_{n+1},\end{array}$$where$$\begin{array}{c}\text{(16)}& W_{n+1}=\begin{array}{c}{w}_{n+p_1}\\ w_{n+p_2}\\ w_{n+p_3}\end{array},W_n=\begin{array}{c}{w}_n\\ w_{n-p_1}\\ w_{n-p_2}\end{array},F_{n+1}=\begin{array}{c}{f}_{n+p_1}\\ f_{n+p_2}\\ f_{n+p_3}\end{array},F_n=\begin{array}{c}{f}_n\\ f_{n-p_1}\\ f_{n-p_2}\end{array},G_{n+1}=\begin{array}{c}{g}_{n+p_1}\\ g_{n+p_2}\\ g_{n+p_3}\end{array},E_1=\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array},E_0=\begin{array}{ccc}1& 0& 0\\ 1& 0& 0\\ 1& 0& 0\end{array},\\ A_1=\begin{array}{ccc}{\alpha }_{1,1}& {\alpha }_{1,2}& {\alpha }_{1,3}\\ {\alpha }_{2,1}& {\alpha }_{2,2}& {\alpha }_{2,3}\\ {\alpha }_{3,1}& {\alpha }_{3,2}& {\alpha }_{3,3}\end{array},A_0=\begin{array}{ccc}{\alpha }_{1,0}& 0& 0\\ {\alpha }_{2,0}& 0& 0\\ {\alpha }_{3,0}& 0& 0\end{array},B_1=\begin{array}{ccc}0& 0& {\beta }_{1,3}\\ 0& 0& {\beta }_{2,3}\\ 0& 0& {\beta }_{3,3}\end{array}.\end{array}$$

3. Analysis of the Method

In this section, we present the truncation error, order, zero-stability and absolute stability of the OMDHBM.

3.1. Local Truncation Error (LTE)

The LTE for the OMDHBM (14) is expressed as$$\begin{array}{}\text{(17)}& {\mathfrak{L}}_{i}w{t}_n;h=\frac{h^{\mu +3}}{\mu +3!}{p}_i^{\mu +3}-\mu +3{\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_i^{\mu +2}-\mu +3\mu +2{\beta }_{i,3}{p}_3^{\mu +1}{w}^{\mu +3}{t}_n+O{h}^{\mu +4}.\end{array}$$

This result was derived by assuming $wt$ to be a sufficiently differentiable function, applying$$\begin{array}{}\text{(18)}& {\mathfrak{L}}_{i}w{t}_n;h=w{t}_n+h{p}_i-w{t}_n-h{\displaystyle \sum }_{j=0}^3{\alpha }_{ij}{w}^{\prime }{t}_n+h{p}_j-h^2{\displaystyle \sum }_{j=1}^3{\beta }_{ij}{w}^{″}{t}_n+h{p}_j,\end{array}$$using Taylor’s series to expand (18) around $t_n$, to obtain$$\begin{array}{}\text{(19)}& \mathfrak{L}w{t}_n;h=U_{0}w{t}_n+U_{1}h{w}^{\prime }{t}_n+U_2{h}^2{w}^{″}{t}_n+U_3{h}^3{w}^{‴}{t}_n+\cdots +U_{\widehat{e}}{h}^{\widehat{e}}{w}^{\widehat{e}}{t}_n,\end{array}$$where $U_1,U_2,U_3,\dots ,U_{\widehat{e}}$ are constants. The OMDHBM have order $\widehat{e}$ if$$\begin{array}{}\text{(20)}& U_0=U_1=\cdots =U_{\widehat{e}}=0,U_{\widehat{e}+1}\ne 0.\end{array}$$

According to [57, 65], expanding equation (19) using the Taylor series, we obtain$$\begin{array}{c}\text{(21)}& {\mathfrak{L}}_{i}w{t}_n;h={\displaystyle \sum }_{\sigma =1}^K\frac{h^{\sigma }{p}_i^{\sigma }}{\sigma !}{w}^{\sigma }{t}_n-{\displaystyle \sum }_{\sigma =1}^K\frac{\sigma h^{\sigma }}{\sigma !}{\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_j^{\sigma -1}{w}^{\sigma }{t}_n-{\displaystyle \sum }_{\sigma =1}^K\frac{\sigma \sigma -1h^{\sigma }}{\sigma !}{\beta }_{i,3}{p}_3^{\sigma -2}{w}^{\sigma }{t}_n+O{h}^{K+1},\end{array}$$

Expanding equation (21), we obtain$$\begin{array}{c}\text{(22)}& {\mathfrak{L}}_{i}w{t}_n;h={\displaystyle \sum }_{\sigma =1}^{\mu +1}\frac{h^{\sigma }{p}_i^{\sigma }}{\sigma !}{w}^{\sigma }{t}_n-{\displaystyle \sum }_{\sigma =1}^{\mu +1}\frac{\sigma h^{\sigma }}{\sigma !}{\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_j^{\sigma -1}{w}^{\sigma }{t}_n-{\displaystyle \sum }_{\sigma =1}^{\mu +1}\frac{\sigma \sigma -1h^{\sigma }}{\sigma !}{\beta }_{i,3}{p}_3^{\sigma -2}{w}^{\sigma }{t}_n+O{h}^{K+1}+{\displaystyle \sum }_{\sigma =\mu +2}^K\frac{h^{\sigma }{p}_i^{\sigma }}{\sigma !}{w}^{\sigma }{t}_n\\ -{\displaystyle \sum }_{\sigma =\mu +2}^K\frac{\sigma h^{\sigma }}{\sigma !}{\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_j^{\sigma -1}{w}^{\sigma }{t}_n-{\displaystyle \sum }_{\sigma =\mu +2}^K\frac{\sigma \sigma -1h^{\sigma }}{\sigma !}{\beta }_{i,3}{p}_3^{\sigma -2}{w}^{\sigma }{t}_n+O{h}^{K+1},\end{array}$$where $\mu =3$ and $K\ge \mu +3$ is a positive integer. This can be expressed as$$\begin{array}{c}\text{(23)}& {\mathfrak{L}}_{i}w{t}_n;h={\displaystyle \sum }_{\sigma =2}^{\mu +1}\frac{h^{\sigma }}{\sigma !}{p}_i^{\sigma }-\sigma {\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_j^{\sigma -1}-\sigma \sigma -1{\beta }_{i,3}{p}_3^{\sigma -2}{w}^{\sigma }{t}_n+{\displaystyle \sum }_{\sigma =\mu +3}^K\frac{h^{\sigma }}{\sigma !}{p}_i^{\sigma }-\sigma {\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_j^{\sigma -1}-\sigma \sigma -1{\beta }_{i,3}{p}_3^{\sigma -2}{w}^{\sigma }{t}_n+O{h}^{K+1}.\end{array}$$

According to [57, 65], it is noted that$$\begin{array}{}\text{(24)}& {\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_j^{\sigma -1}+\sigma -1{\beta }_{i,3}{p}_3^{\sigma -2}=\frac{p_i^{\sigma }}{\sigma },\sigma =2,\dots ,\mu +1.\end{array}$$

Substituting (24) into (23), we obtain$$\begin{array}{}\text{(25)}& {\mathfrak{L}}_{i}w{t}_n;h=\frac{h^{\mu +3}}{\mu +3!}{p}_i^{\mu +3}-\mu +3{\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_i^{\mu +2}-\mu +3\mu +2{\beta }_{i,3}{p}_3^{\mu +1}{w}^{\mu +3}{t}_n+O{h}^{\mu +4}.\end{array}$$

Using equation (25), the error constant is$$\begin{array}{}\text{(26)}& U_{i,\mu +3}=\frac{p_i^{\mu +3}}{\mu +3!}-\frac{1}{\mu +2!}{\displaystyle \sum }_{j=0}^4{\alpha }_{ij}{p}_i^{\mu +2}-\frac{1}{\mu +1}{\beta }_{i3}{p}_i^{\mu +1}.\end{array}$$

To determine a suitable value for $c$, we optimize the LTE in equation (13). This decision is driven by the fact that at the end of the interval [$t_n,t_{n+1}$], only $w_{n+1}$ is required to proceed with computations in the next block. This approach effectively increases the order of accuracy in each equation. The value of $c$ in the LTE corresponding to $w_{n+1}$ is$$\begin{array}{c}\text{(27)}& {\mathfrak{L}}_{i}w{t}_n;h=-\frac{h^{6}5c^2-5c+1w^6{t}_n}{7200}-\frac{h^{7}119c^2-119c+24w^7{t}_n}{302400}-\frac{h^{8}10c^4-20c^3+86c^2-76r+15w^8{t}_n}{604800}\\ -\frac{h^{9}33c^4-66c^3+138c^2-105c+20w^9{t}_n}{3628800}+O{h}^{10},\end{array}$$and is the root of $-1+5c-5c^2$, which can be calculated by equating to zero. We obtained$$\begin{array}{}\text{(28)}& c=\frac{1}{10}5-\sqrt{5}.\end{array}$$

Substituting (28) in (27), the optimized truncation error (27) becomes$$\begin{array}{}\text{(29)}& -\frac{h^7{w}^7{t}_n}{1512000}+O{h}^8.\end{array}$$

The truncation errors for the method (11)–(13) are given below:$$\begin{array}{}\text{(30)}& \frac{5+\sqrt{5}{h}^6{w}^6{t}_n}{450000\sqrt{5}-1}+O{h}^7,\frac{2\sqrt{5}-5h^6{w}^6{t}_n}{225000\sqrt{5}-1}+O{h}^7,-\frac{h^7{w}^7{t}_n}{1512000}+O{h}^8.\end{array}$$

The error constant is$$\begin{array}{}\text{(31)}& U_6={\frac{5+\sqrt{5}}{450000\sqrt{5}-1},\frac{2\sqrt{5}-5}{225000\sqrt{5}-1},0}^T,\end{array}$$which indicates that OMDHBM has an order of $\widehat{e}=5$.

3.2. Zero Stability

The zero stability of the OMDHBM can be determined by considering the matrix (15) as $h⟶0$. In this limit, the matrix (15) reduces to$$\begin{array}{}\text{(32)}& E_1{W}_{n+1}-E_0{W}_n=0.\end{array}$$

The characteristics polynomial is defined as$$\begin{array}{}\text{(33)}& \rho \lambda =\det E_1\lambda -E_0={\lambda }^2\lambda -1.\end{array}$$

By equating $\rho \lambda =0$, the roots of the polynomials are $\lambda =0$ and $\lambda =1$. Thus, following [66], the OMDHBM is zero stable. Since the order of the OMDHBM is greater than 1 ($\widehat{e}=5>1$), the method is consistent. When the method is zero stable and consistent, this leads to the conclusion that the method converges.

3.3. Proof of Absolute Stability

A stable numerical method is affirmed when errors introduced in a single time step do not amplify in subsequent time steps. According to [57, 67], a region of absolute stability for a method is described by$$\begin{array}{}\text{(34)}& \mathfrak{R}\vartheta =\vartheta \mathbb{\in }\mathbb{C}:\mathcal{Q}\vartheta <1.\end{array}$$

A method is considered A-stable if the absolute stability region contains the entire left side of the plane $\mathfrak{R}\vartheta <0$ and L-stable, if ${\lim }_{\vartheta ⟶\infty }\mathfrak{R}\vartheta =0$. The test equations$$\begin{array}{}\text{(35)}& w^{\prime }=\lambda w,w^{″}={\lambda }^{2}w,\end{array}$$are applied to the matrix equation (15) to give$$\begin{array}{}\text{(36)}& W_{n+1}=\mathcal{Q}\vartheta W_n,\vartheta ={\lambda }^2{h}^2,\end{array}$$where $\mathcal{Q}\vartheta$ is a matrix given by$$\begin{array}{}\text{(37)}& \mathcal{Q}\vartheta ={E_1-\vartheta A_1-{\vartheta }^2{B}_1}^{-1}{E}_0+\vartheta A_0.\end{array}$$

By computing the eigenvalues of the matrix $\mathcal{Q}\vartheta$, a stability function $\mathfrak{R}\vartheta$ can be deduced as$$\begin{array}{}\text{(38)}& \mathfrak{R}\vartheta =-\frac{\sqrt{5}+1{\vartheta }^4-4\sqrt{5}{\vartheta }^3-60\sqrt{5}+60{\vartheta }^2+600\sqrt{5}+600\vartheta -6000}{10{\vartheta }^4-140{\vartheta }^3+960{\vartheta }^2-3600\vartheta +6000}.\end{array}$$

As shown in Figure 1, the stability region for the OMDHBM is the shaded area, which encompasses the entire left side of the plane (A-stable), while the unshaded area is the unstable region. It is evident that the method is not L-stable since ${\lim }_{\vartheta ⟶\infty }\mathfrak{R}\vartheta \ne 0$. A-stable method is more suitable for a wider range of problems compared to methods with limited regions of absolute stability. It ensures stability regardless of the step size.

[figure(s) omitted; refer to PDF]

4. Implementation

In this section, we discuss the extension of OMDHBM to solve parabolic PDEs (3). Equation (3) is solved using the OMDHBM with respect to $t$ and the SCM in $x$ after transforming $x\in a,b$ to $\widehat{x}\in -1,1$ using a linear transformation of$$\begin{array}{}\text{(39)}& x=\frac{\widehat{x}}{2}b-a+\frac{1}{2}a+b,\end{array}$$as suggested in [68]. This transformation ensures that the domain of $x$ is mapped onto the interval $-1,1$. The spectral approach replaces the continuous derivatives $w\widehat{x},t$ with discrete approximations at specific grid points, represented as a matrix vector product. The discretized derivatives are expressed as$$\begin{array}{}\text{(40)}& w^{\prime }\widehat{x},t⟶{\mathrm{W}}^{\prime }=\text{DW},w^{″}\widehat{x},t⟶{\mathrm{W}}^{″}={\mathrm{D}}^2\mathrm{W},\end{array}$$where $\mathrm{W}={w{\widehat{x}}_0,t,w{\widehat{x}}_1,t,\dots ,w{\widehat{x}}_{N_x},t}^T$ and $T$ represent matrix transpose. The matrix $D=2D/b-a$, where $D$ is the Chebyshev differentiation matrix [68]. The collocation points used are the Gauss-Lobatto collocation points defined on the interval $-1,1$ as$$\begin{array}{}\text{(41)}& {\widehat{x}}_j=\cos \frac{{\pi }_j}{N_x},j=0,1,\dots ,N_x.\end{array}$$

The application of the spectral method is best described for linear PDEs which can be expressed in general form$$\begin{array}{}\text{(42)}& fx,t,w,w^{\prime },w^{″},\dots ,w^k=\alpha x,t+{\displaystyle \sum }_{\sigma =0}^k{\gamma }_{\sigma }x,t{w}^{\sigma }x,t,\end{array}$$where $\alpha x,t$ and ${\gamma }_{\sigma }x,t$ are known functions. For nonlinear PDEs, the equations are linearized before applying the SCM. Utilizing the SCM with $N_x+1$ collocation points on equation (3) gives equation (5) of the form$$\begin{array}{}\text{(43)}& \dot{w}=ft,w,w0=w_0,\end{array}$$with BCs given by$$\begin{array}{}\text{(44)}& w{\widehat{x}}_N,t={\aleph }_{a}t,w{\widehat{x}}_0,t={\aleph }_{b}t.\end{array}$$

Utilizing the OMDHBM on equation (43) gives$$\begin{array}{}\text{(45)}& w_{n+p_i}=w_n+h{\gamma }_i{f}_n+h{\displaystyle \sum }_{j=1}^3{\alpha }_{ij}{f}_{n+p_j}+h^2{\displaystyle \sum }_{j=3}^1{\beta }_{ij}{g}_{n+p_j},i=1,2,3,\end{array}$$where$$\begin{array}{}\text{(46)}& f_{n+p_j}=f{t}_{n+p_j},w_{n+p_j},f_n=f{t}_n,w_n,g_{n+p_j}=g{t}_{n+p_j},w_{n+p_j},f_{n+p_j}.\end{array}$$

Equation (45) can be defined as$$\begin{array}{c}\text{(47)}& w_{n+p_i}=w_n+h{\gamma }_i{f}_n+h{\displaystyle \sum }_{j=1}^3{\alpha }_{ij}{\mathrm{\Phi }}_{n+p_j}{w}_{n+p_j}+{\mathrm{\Psi }}_{n+p_j}+h^2{\displaystyle \sum }_{j=3}^1{\beta }_{ij}{K}_{n+p_j}+{\mathrm{\Omega }}_{n+p_j}{w}_{n+p_j}+{\mathrm{\Theta }}_{n+p_j}{\mathrm{\Phi }}_{n+p_j}{w}_{n+p_j}+{\mathrm{\Psi }}_{n+p_j},\end{array}$$where $i=1,2,3.$ In matrix form, equation (47) becomes$$\begin{array}{c}\text{(48)}& \begin{array}{ccc}{E}_{11}& E_{12}& E_{13}\\ E_{21}& E_{22}& E_{23}\\ E_{31}& E_{32}& E_{33}\end{array}\begin{array}{c}{w}_{n+p_1}\\ w_{n+p_2}\\ w_{n+p_3}\end{array}=\begin{array}{c}{w}_n\\ w_n\\ w_n\end{array}+\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{21}& A_{22}& A_{23}\\ A_{31}& A_{32}& A_{33}\end{array}\begin{array}{c}{f}_n\\ f_n\\ f_n\end{array}+\begin{array}{ccc}{C}_{11}& C_{12}& C_{13}\\ C_{21}& C_{22}& C_{23}\\ C_{31}& C_{32}& C_{33}\end{array}\begin{array}{c}{\mathrm{\Psi }}_{n+p_1}\\ {\mathrm{\Psi }}_{n+p_2}\\ {\mathrm{\Psi }}_{n+p_3}\end{array}\\ +\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{21}& B_{22}& B_{23}\\ B_{31}& B_{32}& B_{33}\end{array}\begin{array}{c}{K}_{n+p_1}\\ K_{n+p_2}\\ K_{n+p_3}\end{array}\\ +\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{21}& B_{22}& B_{23}\\ B_{31}& B_{32}& B_{33}\end{array}\begin{array}{ccc}{\mathrm{\Theta }}_{n+p_1}& 0& 0\\ 0& {\mathrm{\Theta }}_{n+p_2}& 0\\ 0& 0& {\mathrm{\Theta }}_{n+p_3}\end{array}\begin{array}{c}{\mathrm{\Psi }}_{n+p_1}\\ {\mathrm{\Psi }}_{n+p_2}\\ {\mathrm{\Psi }}_{n+p_3}\end{array},\end{array}$$where$$\begin{array}{c}\text{(49)}& {\mathrm{\Phi }}_{n+p_j}=H_{2,n+p_j}t{D}^2+H_{1,n+p_j}tD+H_{0,n+p_j}t,{\mathrm{\Omega }}_{n+p_j}=Q_{2,n+p_j}t{D}^2+Q_{1,n+p_j}tD+Q_{0,n+p_3}t,\\ {\mathrm{\Theta }}_{n+p_j}=S_{2,n+p_j}t{D}^2+S_{1,n+p_j}tD+S_{0,n+p_j}t,\end{array}$$where ${\mathrm{\Phi }}_{n+p_j}$, ${\mathrm{\Omega }}_{n+p_j}$ and ${\mathrm{\Theta }}_{n+p_j}$ are diagonal matrices.$$\begin{array}{c}\text{(50)}& E_{ij}=\begin{array}{cc}-h{\alpha }_{i,j}{\mathrm{\Phi }}_{n+pj}-h^2{\beta }_{i,j}{\mathrm{\Omega }}_{n+p_j}-h^2{\beta }_{i,j}{\mathrm{\Theta }}_{n+p_j}{\mathrm{\Phi }}_{n+p_j},& i\ne j,\\ I-h{\alpha }_{i,j}{\mathrm{\Phi }}_{n+pj}-h^2{\beta }_{i,j}{\mathrm{\Omega }}_{n+p_j}-h^2{\beta }_{i,j}{\mathrm{\Theta }}_{n+p_j}{\mathrm{\Phi }}_{n+p_j},& i=j,\end{array}{A}_{ij}=h{\gamma }_{i,j}I,C_{ij}=h{\alpha }_{i,j}I,B_{ij}=h^2{\beta }_{i,j}I,\end{array}$$where $I$ is an identity matrix of order $N_x+1$, $N_x$ is the number of intervals and$$\begin{array}{c}\text{(51)}& {\alpha }_{ij}=\begin{array}{ccc}\frac{1}{600}155-11\sqrt{5}& \frac{1}{600}95-71\sqrt{5}& \frac{97\sqrt{5}-25}{3000}\\ \frac{1}{600}71\sqrt{5}+95& \frac{1}{600}11\sqrt{5}+155& \frac{-97\sqrt{5}-25}{3000}\\ \frac{5}{12}& \frac{5}{12}& \frac{1}{12}\\ \hspace{1em}& \hspace{1em}& \hspace{1em}\end{array},{\gamma }_{ij}=\begin{array}{ccc}\frac{13\sqrt{5}+275}{3000}& 0& 0\\ 0& \frac{275-13\sqrt{5}}{3000}& 0\\ 0& 0& \frac{1}{12}\end{array},{\beta }_{ij}=\begin{array}{ccc}0& 0& -\frac{1}{50\sqrt{5}}\\ 0& 0& \frac{1}{50\sqrt{5}}\\ 0& 0& 0\\ \hspace{1em}& \hspace{1em}& \hspace{1em}\end{array}.\end{array}$$

The given expression in (48) is expressed compactly as$$\begin{array}{}\text{(52)}& E_1{W}_{n+1}=W_n+A_0{F}_n+C_1{\mathrm{\Psi }}_{n+1}+B_1{K}_{n+1}+B_1{\mathrm{\Theta }}_{n+1}{\mathrm{\Psi }}_{n+1}.\end{array}$$

This can also be expressed in the form$$\begin{array}{}\text{(53)}& E_1{W}_{n+1}=R_{n+1},\end{array}$$where $R_{n+1}=W_n+A_0{F}_n+C_1{\mathrm{\Psi }}_{n+1}+B_1{K}_{n+1}+B_1{\mathrm{\Theta }}_{n+1}{\mathrm{\Psi }}_{n+1}$.

The BCs for the PDE are imposed on the first and last row of each $E_1$ and the corresponding row entry of $R_{n+1}$ in (53) as follows:$$\begin{array}{}\text{(54)}& \begin{array}{c}\begin{array}{ccc}\begin{array}{cc}1& 0\end{array}& \cdots & \begin{array}{cc}0& 0\end{array}\end{array}\\ E_1\\ \begin{array}{ccc}\begin{array}{cc}0& 0\end{array}& \dots & \begin{array}{cc}0& 1\end{array}\end{array}\end{array}\begin{array}{c}{w}_{n+p_1}\\ w_{n+p_2}\\ w_{n+p_3}\end{array}=\begin{array}{c}wb,t\\ R\\ wa,t\end{array}.\end{array}$$

The next section discusses the applications of the SOMDHBM on the F-N equations (1) and (2).

5. Applications in F-N Equations

Prior to applying the method, the F-N equations are linearized using a linear partition technique. The core concept behind this linearization approach is to evaluate all linear terms at the present iteration level $r+1$ and all nonlinear terms at the previous iteration level $r$. The difference between the iterations is assumed to be very small. The details of the linear partition technique are provided in [44, 69]. The linearized form of the F-N equation (1) is$$\begin{array}{}\text{(55)}& {\dot{w}}_{r+1}=w_{r+1}^{″}-\delta w_{r+1}+w_r^{2}1+\delta -w_r.\end{array}$$

The second-derivative with respect to $t$ is$$\begin{array}{c}\text{(56)}& {\ddot{w}}_{r+1}=-\delta w_{r+1}^{″}+{\delta }^2{w}_{r+1}+w_{r^{″}}{w}_{r^{″}}+2+\delta w_r+\delta -2w_r^2-w_r^3-3\delta +3{\delta }^2{w}_r^2+2{\delta }^2+8\delta +2w_r^3-5+5\delta w_r^4+3w_r^5.\end{array}$$

The SOMDHBM is applied using the parameters from equation (55)$$\begin{array}{}\text{(57)}& {\mathbf{\Phi }}_{n+p_j}={\mathit{D}}^2-\delta \mathit{I},{\mathbf{\Psi }}_{n+p_j}=w_r^{2}1+\delta -w_r,\end{array}$$and from equation (56)$$\begin{array}{c}\text{(58)}& {\mathit{\Theta }}_{n+p_j}=0,{\mathit{\Omega }}_{n+p_j}=-\delta {\mathit{D}}^2+{\delta }^2\mathit{I},{\mathit{K}}_{n+p_j}={\mathit{D}}^2{w}_r+2+\delta w_r+\delta -2w_r^2-w_r^3\\ {\mathit{D}}^2{w}_r-3\delta +3{\delta }^2{w}_r^2+2{\delta }^2+8\delta +2w_r^3-5+5\delta w_r^4+3w_r^5,\end{array}$$where $\mathit{I}$ is $N_x+1$ × $N_x+1$ identity matrix and $N_x$ is the number of intervals. The step-by-step implementation of the SOMDHBM for solving the F-N equations is presented in Algorithm 2. All numerical simulations were performed using MATLAB 2016a on a 64 bit Windows system with an Intel(R) Core i5 processor (2 physical cores, 4 logical cores, 2.50 GHz), 8GB RAM. The computations were executed using 2 logical cores, as assigned by MATLAB due to hyper-threading.

Algorithm 2: Algorithmic implementation of the SOMDHBM.

  Input: Parameter $\mathit{\delta }$, step size $\mathit{h}$, number of spatial points ${\mathit{N}}_{\mathit{x}}$, number of iterations $\mathit{R}$

1. Generate Chebyshev collocation points and compute $\mathit{D}$, ${\mathit{D}}_{\mathbf{1}}$, ${\mathit{D}}_{\mathbf{2}}$ over spatial domain $\mathit{a},\mathit{b}$;

2. Define the time grid ${\mathit{N}}_{\mathit{t}}$ with step size $\mathit{h}$ over an interval $\mathbf{0},\mathit{T}$;

3. Define initial condition $\mathit{w}\mathit{x},\mathbf{0}$;

4. Define boundary functions $\mathit{w}\mathit{a},\mathit{t}$, $\mathit{w}\mathit{b},\mathit{t}$, and exact solution $\mathit{w}\mathit{x},\mathit{t}$;

5. Define first-derivative parameters: $\mathrm{\Phi }$, $\mathrm{\Psi }$;

6. Define second-derivative parameters: $\mathrm{\Theta }$, $\mathrm{\Omega }$, $\mathit{K}$;

7. For$\mathit{n}=\mathbf{1}$ to ${\mathit{N}}_{\mathit{t}}$:

   a. Initialize guesses: ${\mathit{W}}_{\mathit{r}}:,\mathit{n}+\mathit{i}\leftarrow \mathit{W}:,\mathit{n}$ for $\mathit{i}=\mathbf{0}$ to $\mathbf{3}$;

   b. For$\mathit{r}=\mathbf{1}$ to $\mathit{R}$:

    i. For$\mathit{i}=\mathbf{0}$ to $\mathbf{3}$:

     - Evaluate $\mathrm{\Psi }$, $\mathrm{\Phi }$, $\mathrm{\Theta }$, $\mathrm{\Omega }$, $\mathit{K}$ at ${\mathit{t}}_{\mathit{n}}+\mathit{h}\mathbf{·}\text{points}\mathit{i}+\mathbf{1}$;

    ii. For$\mathit{i}=\mathbf{1}$ to $\mathbf{3}$:

     - For$\mathit{j}=\mathbf{1}$ to $\mathbf{3}$:

     ∗ Assemble coefficient matrices ${\mathit{E}}_{\mathit{i},\mathit{j}}$, ${\mathit{A}}_{\mathit{i},\mathit{j}}$, ${\mathit{C}}_{\mathit{i},\mathit{j}}$, ${\mathit{B}}_{\mathit{i},\mathit{j}}$;

     ∗ Impose boundary conditions on rows 1 and ${\mathit{N}}_{\mathit{x}}+\mathbf{1}$ of ${\mathit{E}}_{\mathit{i},\mathit{j}}$;

     - Adjust diagonals in ${\mathit{E}}_{\mathit{i},\mathit{i}}$ for boundary rows;

    iii. Flatten block matrices into matrix system: ${\mathit{E}}_{\mathbf{1}}\mathbf{·}\mathit{W}={\mathit{W}}_{\mathit{n}}+{\mathit{A}}_{\mathbf{0}}\mathbf{·}{\mathit{F}}_{\mathit{n}}+{\mathit{C}}_{\mathbf{1}}\mathbf{·}\mathbf{\Psi }+{\mathit{B}}_{\mathbf{1}}\mathbf{·}\mathit{K}+{\mathit{B}}_{\mathbf{1}}\mathbf{·}\mathbf{\Theta }\mathbf{·}\mathbf{\Psi }$;

    iv. Apply boundary values in RHS vector;

    v. Solve the linear system for ${\mathit{W}}_{\mathit{n}+\mathit{i},\mathit{r}}$;

    vi. Compute error: $\parallel {\mathit{W}}_{\mathit{n}+\mathbf{3},\mathit{r}}-\text{exact}\mathit{x},{\mathit{t}}_{\mathit{n}}+\mathit{h}\mathbf{·}\text{points}\mathbf{4}{\parallel }_{\mathbf{\infty }}$;

   c. Update solution: $\mathit{W}:,\mathit{n}+\mathbf{1}\leftarrow {\mathit{W}}_{\mathit{r}}:,\mathit{n}+\mathbf{3}$;

  Output: Final numerical solution $\mathit{W}$

Example 1.

Consider the F-N equation (1) in [62], subject to the BCs$$\begin{array}{c}\text{(59)}& wa,t=\frac{\delta }{1+e^{-2A_{1}a-A_{2}t}},wb,t=\frac{\delta }{1+e^{-2A_{1}b-A_{2}t}},t\ge 0,\end{array}$$and the IC$$\begin{array}{}\text{(60)}& wx,0=\frac{\delta }{1+e^{-2A_{1}x}},x\in 0,1.\end{array}$$

The exact solution of (1) is$$\begin{array}{}\text{(61)}& wx,t=\frac{\delta }{1+e^{-2A_{1}x-A_{2}t}},x\in 0,1,t\ge 0,\end{array}$$where$$\begin{array}{}\text{(62)}& A_1=\frac{\sqrt{2}}{4}\delta ,A_2=\frac{4-2\delta }{\delta }{A}_1.\end{array}$$

The BCs are imposed on the first and last rows of equation (53) in the form of equation (54). For different values of $\delta$, the SOMDHBM for the current problem is evaluated and compared to the calculated absolute error for accuracy and efficiency. The absolute error is defined as$$\begin{array}{}\text{(63)}& A{E}_n=w_{n+p_i}-W_{n+p_i},\end{array}$$where $w_{n+p_i}$ and $W_{n+p_i}$ represent the exact and numerical solutions, respectively.

When $\delta =0.01,N_x=50,h=0.01$ and fifteen iterations, Table 1 shows a comparison of the absolute error ($A{E}_n$) in the computed solutions for the F-N equation (1) obtained using ADM [11], VIM [12, 13], DTM [15], NSFD [62] and SOMDHBM. The results demonstrate that SOMDHBM consistently achieves the lowest $A{E}_n$ error across all test cases, with errors significantly lower than those obtained using other methods. For instance, at $x=0.1$ and $t=1.0$, the SOMDHBM achieves an $A{E}_n$ error in the order of ${10}^{-13}$, whereas NSFD produces errors around ${10}^{-12}$ and ADM, VIM and DTM yield errors in the range of ${10}^{-6}$ to ${10}^{-5}$. This confirms the high accuracy of SOMDHBM, as it outperforms all compared methods. The $\text{NFE}val$ column in table indicates the number of function evaluations required by SOMDHBM. The values show that SOMDHBM maintains its high accuracy while using a relatively small number of function evaluations. Based on the observed numerical results, the SOMDHBM is highly competitive with existing methods and demonstrates efficiency in finding solutions for the F-N equation (1).

When $\delta =0.001,N_x=50,h=0.01$ and fifteen iterations, Table 2 presents the $A{E}_n$ error for different values of $x$ and $t$ in the F-N (1), comparing the performance of the SOMDHBM and NSFD [62]. The data reveal that SOMDHBM consistently produces more accurate results than NSFD, achieving significantly lower $A{E}_n$ error. At all tested spatial points ($x=0.2,0.5,0.8$), the $A{E}_n$ error of SOMDHBM remain within the order of ${10}^{-15}$, while NSFD exhibits errors in the order of ${10}^{-14}$, which are approximately ten times larger. This trend is observed across different time steps, confirming that SOMDHBM maintains its precision over time. The $\text{NFE}val$ for SOMDHBM increases as $t$ progresses, from 1800 at $t=0.3$ to 5400 at $t=0.9$. This suggests that while computational effort increases over time, SOMDHBM retains its superior accuracy. These results show that the SOMDHBM not only improves accuracy compared to NSFD but also provides a stable and efficient approach for solving the F-N (1).

The numerical results demonstrate that the SOMDHBM accurately captures the propagation of traveling wavefronts, which represent electrical signals, across a spatial domain over time. The extremely low error magnitudes (on the order of ${10}^{-15}$) indicate that SOMDHBM preserves sharp wave profiles, essential for realistic simulations of nerve activity. The ability of the method to handle small $\delta$ values (which correspond to near-threshold excitability) shows its robustness in modeling sensitive nerve dynamics.

The numerical and exact solutions of (1) for $\delta =0.001$ with $h=0.01$ are compared and presented in Figure 2. The graph shows that the numerical solution aligns closely with the exact solution throughout the entire domain, highlighting the effectiveness of the SOMDHBM. Figure 3 shows the exact solution alongside the numerical solution. The convergence and error profile of the numerical solution for SOMDHBM when $\delta =0.01,0.001$, $h=0.01$ and $N_x=50$ are illustrated in Figure 4. As the $\delta$ decreases, the SOMDHBM converges faster.

Table 1

${AE}_n$ of the solution for F-N equation compared with ADM, VIM, DTM and NSFD, with $\delta =0.01$$N_x=50$ and $h=0.01$.

$x$	$t$	SOMDHBM	$\text{NFE}val$	NSFD [62]	ADM [11]	VIM [12, 13]	DTM [15]
0.1	0.05	$3.4425E-13$	300	$5.0915E-12$	$2.4853E-06$	$2.4880E-06$	$2.4875E-06$
	0.1	$1.3896E-13$	600	$7.4482E-12$	$4.9728E-06$	$4.9750E-06$	$4.9750E-06$
	1.0	$2.0949E-13$	6000	$1.1137E-12$	$4.9744E-05$	$4.9750E-05$	$4.9750E-05$
0.5	0.05	$6.9701E-15$	300	$1.1430E-11$	$2.4763E-06$	$2.4875E-06$	$2.4875E-06$
	0.1	$6.5494E-15$	600	$1.8999E-11$	$4.9638E-06$	$4.9750E-06$	$4.9750E-06$
	1.0	$6.2910E-15$	6000	$3.0936E-11$	$4.9738E-05$	$4.9749E-05$	$4.9749E-05$
0.9	0.05	$3.1650E-15$	300	$5.0915E-12$	$2.4674E-06$	$2.4875E-06$	$2.4874E-06$
	0.1	$3.3081E-15$	600	$7.4481E-12$	$4.9549E-06$	$4.9750E-06$	$4.9750E-06$
	1.0	$3.4564E-15$	6000	$7.8174E-12$	$4.9730E-05$	$4.9789E-05$	$4.9749E-05$

Table 2

${AE}_n$ of the solution for F-N equation compared with NSFD, with $N_x=50$, $\delta =0.001$ and $h=0.01$.

$t$	SOMDHBM				NSFD [62]
	$x=0.2$	$x=0.5$	$x=0.8$	$\text{NFE}val$	$x=0.2$	$x=0.5$	$x=0.8$
0.3	$4.6113E-15$	$6.0347E-16$	$4.2761E-16$	1800	$1.8987E-14$	$1.9838E-14$	$1.9974E-14$
0.5	$4.4390E-15$	$5.8384E-16$	$4.2696E-16$	3000	$2.9530E-14$	$3.0978E-14$	$3.1209E-14$
0.9	$4.2354E-15$	$5.2259E-16$	$4.2577E-16$	5400	$1.8987E-14$	$1.9838E-14$	$1.9974E-14$

[figure(s) omitted; refer to PDF]

Example 2.

Consider the F-N (1) in [49], subject to the BCs$$\begin{array}{c}\text{(64)}& wa,t=\frac{1}{2}+\frac{1}{2}\tanh \frac{1}{2\sqrt{2}}a-\frac{2\delta -1}{\sqrt{2}}t,wb,t=\frac{1}{2}+\frac{1}{2}\tanh \frac{1}{2\sqrt{2}}b-\frac{2\delta -1}{\sqrt{2}}t,\end{array}$$and the IC$$\begin{array}{}\text{(65)}& wx,0=\frac{1}{2}+\frac{1}{2}\tanh \frac{x}{2\sqrt{2}},x\in a,b.\end{array}$$

The exact solution is$$\begin{array}{}\text{(66)}& wx,t=\frac{1}{2}+\frac{1}{2}\tanh \frac{1}{2\sqrt{2}}x-\frac{2\delta -1}{\sqrt{2}}t.\end{array}$$

By imposing the BCs from Example 2 on (53) in the form of (54), we analyze the numerical accuracy and efficiency of different methods. With $a=-10,b=10,\delta =0.75,N_x=50,h=0.02$ over fifteen iterations, Table 3 presents the maximum absolute error (${\text{ME}}_{\infty }$), error norm (${\text{EN}}_2$), computation time and number of function evaluations ($\text{NFE}val$) for different values of $t$ in the F-N equation (1). The table compares the performance of SOMDHBM, PDQM [61] and CFDBM [49].

The numerical results show that SOMDHBM consistently achieves the lowest ${\text{ME}}_{\infty }$ error ranging from ${10}^{-9}$ to ${10}^{-8}$, while PDQM exhibits much higher errors, reaching values as large as ${10}^{-3}$ to ${10}^{-4}$. Similarly, the ${\text{EN}}_2$ norm for SOMDHBM remains significantly smaller than that of PDQM, confirming the effectiveness of the SOMDHBM. Compared to CFDBM, SOMDHBM achieves at least one order of magnitude lower error across all time steps. The SOMDHBM demonstrates higher computational efficiency by achieving lower computation times while maintaining superior accuracy.

The CPU time for SOMDHBM remains below 0.2 s even for large values of $t$, whereas PDQM requires significantly longer processing time, reaching 0.55 s at $t=5.0$. This indicates that SOMDHBM not only produces more accurate results but also operates with significantly lower computation time, making it suitable for long-term and large-scale simulations in applied fields like neuroscience and chemical reaction-diffusion systems. The numerical findings confirm that the SOMDHBM outperforms PDQM and CFDBM in both accuracy and efficiency, making it a highly effective method for solving the F-N equation (1). The maximum absolute error is calculated using the following equation:$$\begin{array}{}\text{(67)}& M{E}_{\infty }=\underset{1\le i\le M}{\max }{w}_{n+p_i}-W_{n+p_i},\end{array}$$and the error norm is$$\begin{array}{}\text{(68)}& E{N}_2={{\displaystyle \sum }_{i=1}^M{w_{n+p_i}-W_{n+p_i}}^2}^{1/2}.\end{array}$$

Figure 5 illustrates the comparison between the numerical and exact solutions of the F-N equation (1) for $\delta =0.75$ with a step-size of $h=0.02$. Figure 6 presents a proximity of the exact solution against the numerical solution. The convergence and corresponding error profile of the SOMDHBM are displayed in Figure 7, indicating that the numerical method is accurate and that convergence is achieved after the 3rd iteration. The efficiency curves in Figure 8 illustrate the relationship between CPU time and error measures (${\text{ME}}_{\infty }$ and ${\text{EN}}_2$) for the SOMDHBM. The plot suggests that as CPU time increases, the error decreases, indicating the accuracy of the method improves with more computational effort.

Table 3

${ME}_{\infty }$ and ${EN}_2$ of the solutions for Example 2 with $\delta =0.75$, $N_x=50$ and $h=0.02$.

$t$	SOMDHBM				CFDBM [49]		PDQM [61]
	${\text{ME}}_{\infty }$	${\text{EN}}_2$	CPU time	$\text{NFE}val$	${\text{ME}}_{\infty }$	${\text{ME}}_{\infty }$	${\text{EN}}_2$	CPU time
0.2	$1.0905E-09$	$6.8424E-10$	0.02	600	$5.4996E-08$	$4.7416E-05$	$2.3012E-06$	0.10
0.5	$2.8575E-09$	$1.1878E-09$	0.02	1500	$1.0653E-07$	$1.2312E-04$	$5.5695E-06$	0.10
1.0	$6.1097E-09$	$1.8586E-09$	0.04	3000	$1.5822E-07$	$2.6261E-04$	$1.1864E-05$	0.18
1.5	$9.7815E-09$	$2.4373E-09$	0.05	4500	$1.9190E-07$	$4.2096E-04$	$1.9400E-05$	0.24
2.0	$1.3928E-08$	$2.9581E-09$	0.07	6000	$2.1710E-07$	$5.9999E-04$	$2.8162E-05$	0.31
3.0	$2.3928E-08$	$3.9140E-09$	0.11	9000	$2.5517E-07$	$1.0324E-03$	$4.9735E-04$	0.35
5.0	$5.3186E-08$	$5.7975E-09$	0.17	15,000	$3.1252E-07$	$2.3050E-03$	$1.1395E-04$	0.55

[figure(s) omitted; refer to PDF]

Example 3.

Consider the F-N equation (1), as presented in [49], subject to the BCs$$\begin{array}{}\text{(69)}& wa,t=\frac{1}{1+e^{-a-A_{1}t/\sqrt{2}}},wb,t=\frac{1}{1+e^{-b-A_{1}t/\sqrt{2}}},\end{array}$$and the IC$$\begin{array}{}\text{(70)}& wx,0=\frac{1}{1+e^{-x/\sqrt{2}}}.\end{array}$$

The exact solution is$$\begin{array}{}\text{(71)}& wx,t=\frac{1}{1+e^{-A_2/\sqrt{2}}},x\in 0,1,t\ge 0,\end{array}$$where$$\begin{array}{}\text{(72)}& A_1=\sqrt{2}\frac{1}{2}-\delta ,A_2=x+A_{1}t.\end{array}$$

Table 4 displays the absolute error (${\text{AE}}_n$) and computation time for different values of $x$ at $t=0.04$ over one iteration for $N=10$ in equation (1). It compares the performance of the SOMDHBM with CFDBM [49], İnan $\text{et\hspace{0.17em}al}.$ [63] and ExpFDM [64]. The results indicate that SOMDHBM achieves the highest accuracy, with ${\text{AE}}_n$ error ranging from ${10}^{-15}$ to ${10}^{-16}$, while CFDBM produces errors on the order of ${10}^{-9}$ to ${10}^{-8}$. Although CFDBM is more accurate than İnan $\text{et\hspace{0.17em}al}.$ and ExpFDM, SOMDHBM still outperforms CFDBM by at least six orders of magnitude across all spatial points. The SOMDHBM, like CFDBM, converges in just one iteration, whereas İnan $\text{et\hspace{0.17em}al}.$ and ExpFDM require eight iterations, reinforcing the computational efficiency of the SOMDHBM.

Table 5 further supports these findings by presenting the $\text{NFE}val$, ${\text{ME}}_{\infty }$ error and ${\text{EN}}_2$ norm for SOMDHBM. The results show that ${\text{ME}}_{\infty }$ and ${\text{EN}}_2$ remain stable across all $x$, with values between ${10}^{-14}$ and ${10}^{-15}$, confirming that SOMDHBM maintains its accuracy across the domain. Overall, the numerical findings confirm that SOMDHBM surpasses CFDBM, İnan $\text{et\hspace{0.17em}al}.$, and ExpFDM in both accuracy and efficiency, making it an optimal approach for solving the F-N equation (1). Its ability to capture sharp gradients and preserve wave profiles with minimal iterations also makes it ideal for real-time simulations and high-fidelity modeling, enhancing its practical relevance in both theoretical and applied contexts.

Figure 9 compares the numerical solution to the exact solution of the F-N equation (1) when $\delta =0.5,h=0.02$ over fifteen iterations. Figure 10 shows the closeness of the exact solution and the numerical solution. The convergence and corresponding error profile of the SOMDHBM are depicted in Figure 11. The efficiency curves in Figure 12 show the relationship between CPU time, ${\text{ME}}_{\infty }$ and ${\text{EN}}_2$ for the SOMDHBM. The graph shows that the method efficiently reduces computational errors with minimal processing time.

Table 4

${AE}_n$ of the solution for Example 3, with $\delta =0.5,t=0.04$ and $N=10$.

$x$	SOMDHBM	CPU time	CFDBM [49]	İnan $\text{et\hspace{0.17em}al}.$ [63]	ExpFDM [64]
0.2	$4.3299E-15$	0.0014	$1.72E-08$	$2.00E-07$	$3.00E-06$
0.4	$7.7716E-16$	0.0016	$4.63E-09$	$5.00E-07$	$1.00E-05$
0.6	$7.7716E-16$	0.0016	$1.63E-09$	$7.00E-07$	$2.00E-05$
0.8	$1.2212E-15$	0.0016	$7.26E-09$	$6.00E-07$	$4.00E-05$
1	$3.1086E-15$	0.0022	0	0	0
Iterations	1		1	8	8

Table 5

${ME}_{\infty }$ and ${EN}_2$ of the solution for Example 3, with $\delta =0.5,t=0.04$ and $N=10$.

SOMDHBM
$x$	${\text{ME}}_{\infty }$	${\text{EN}}_2$	$\text{NFE}val$
0.2	$4.1078E-14$	$2.6128E-14$	40
0.4	$6.9944E-15$	$3.9893E-15$	40
0.6	$6.2172E-15$	$3.8373E-15$	40
0.8	$9.9920E-15$	$5.3106E-15$	40
1	$3.0420E-14$	$1.4457E-14$	40
Iteration	1

[figure(s) omitted; refer to PDF]

Example 4.

Consider the F-N equation (2), with the exact solutions given by$$\begin{array}{c}\text{(73)}& wx,t=e^{-\epsilon \mathrm{\Delta }t/3}{A}_1{e}^{-cx}+A_2{e}^{cx},x\in -10,10,t\ge 0,\text{(74)}& yx,t=-\frac{1}{6\mathrm{\Delta }{e^{\epsilon \mathrm{\Delta }t/3}}^3{e^{cx}}^3}2{e^{cx}}^6{A_2}^3\mathrm{\Delta }+6{e^{cx}}^4{A}_1{A_2}^2\mathrm{\Delta }-9{e^{cx}}^4{e^{\epsilon \mathrm{\Delta }t/3}}^2{A}_2+6{e^{cx}}^2{A_1}^2{A}_2\mathrm{\Delta }-9{e^{cx}}^2{e^{\epsilon \mathrm{\Delta }t/3}}^2{A}_1+2{A_1}^3\mathrm{\Delta },\end{array}$$where$$\begin{array}{}\text{(75)}& c=\frac{\sqrt{-2{\epsilon \mathrm{\Delta }}^4-6{\epsilon \mathrm{\Delta }}^3+9{\epsilon }^3{\mathrm{\Delta }}^2}\sqrt{6}}{6\epsilon \mathrm{\Delta }\sqrt{\mathrm{\Lambda }}\sqrt{\epsilon \mathrm{\Delta }}},\end{array}$$

as presented in [9, 10].

With parameters $A_1=1,A_2=1,\epsilon =0.3,\mathrm{\Lambda }=1.03,\mathrm{\Delta }=2,N_x=50,$ and $h=0.1$, Table 6 presents the accuracy and efficiency of the SOMDHBM in solving a nonlinear neuronal PDE system. The method achieves very high accuracy, with ${\text{ME}}_{\infty }$ and ${\text{EN}}_2$ remaining extremely small (around ${10}^{-13}$ to ${10}^{-12}$) across time steps. The errors increase slowly and smoothly, indicating good long-term stability. Computationally, the method is efficient, with CPU time rising moderately and the number of function evaluations increasing linearly. Overall, SOMDHBM is shown to be accurate, stable, and reliable for solving complex neurodynamic systems.

Figures 13, 14, and 15 illustrate the comparison between the numerical and exact solutions of the variables $w$ and $y$ for Example 4, using parameters $N_x=50,h=0.1$ and $t=10.$ Figure 13 displays both the numerical and exact solutions of $w$ and $y$ over the spatial domain $x\in -10,10.$, demonstrating a close agreement and validating the accuracy of SOMDHBM. Figures 14 and 15 separately present the numerical and exact solutions for $w$ and $y$, respectively. The near-perfect alignment between the numerical and exact profiles confirms the reliable performance of the method for each component. Figure 16 illustrates the convergence behavior of the numerical solutions $w$ and $y$ over fifteen iterations. After approximately the sixth iteration, the norms of both solutions decrease rapidly and nearly exponentially, confirming strong stability and rapid convergence of the SOMDHBM.

Table 6

${ME}_{\infty }$ and ${EN}_2$ of the SOMDHBM solutions for Example 4 with $N_x=50$ and $h=0.1$.

$t$	$w$		$y$		CPU time	$\text{NFE}val$
	${\text{ME}}_{\infty }$	${\text{EN}}_2$	${\text{ME}}_{\infty }$	${\text{EN}}_2$
1.0	$6.8345E-13$	$9.2908E-13$	$9.2584E-13$	$1.4409E-12$	0.31	750
2.0	$1.1466E-12$	$1.1409E-12$	$2.4992E-12$	$2.6659E-12$	0.40	1500
3.0	$1.2450E-12$	$1.1523E-12$	$3.9191E-12$	$3.3491E-12$	0.29	2250
4.0	$1.5676E-12$	$1.2269E-12$	$4.9390E-12$	$3.6615E-12$	0.45	3000
5.0	$2.2278E-12$	$1.4265E-12$	$5.5462E-12$	$3.7777E-12$	0.46	3750
6.0	$3.1541E-12$	$1.7283E-12$	$5.9411E-12$	$3.8094E-12$	0.56	4500
7.0	$4.3344E-12$	$2.1448E-12$	$6.2205E-12$	$3.8135E-12$	0.66	5250
8.0	$5.7888E-12$	$2.8069E-12$	$6.4206E-12$	$3.8143E-12$	0.88	6000
9.0	$8.0435E-12$	$4.0106E-12$	$6.6980e-12$	$3.8204E-12$	0.89	6750
10.0	$1.1832E-11$	$6.1989E-12$	$7.6159E-12$	$3.8380E-12$	0.89	7500

[figure(s) omitted; refer to PDF]

6. Conclusions

In this study, the efficiency and accuracy of the SOMDHBM were demonstrated for the F-N equation over large space and time intervals. The method utilized the optimized MDHBM in time and the SCM in space to solve F-N equation, a widely studied model in neuroscience and mathematical biology. The SOMDHBM results were compared with those obtained from the ADM, İnan $\text{et\hspace{0.17em}al}.$ [63], CFDBM, ExpFDM, DTM, NSFD, PDQM and VIM scheme. Details of the numerical results were presented and it was evident from the numerical experiments that the SOMDHBM solutions performed better than other schemes in terms of accuracy, computation efficiency and agrees with the exact solutions. Future work may extend this approach to higher-dimensional or stochastic versions of the F-N model, or apply it to other nonlinear PDEs arising in applied sciences.

Funding

This research was funded by the University of KwaZulu-Natal, South Africa.

Acknowledgments

The authors are grateful to the University of Kwazulu-Natal, South Africa, for assistance.

Glossary

Nomenclature

ADMAdomian decomposition method

AEAbsolute error

BCsBoundary conditions

BHMBlock hybrid method

CFDBMCompact finite difference block method

CPU timeComputation time

DQMDifference quadrature method

DTMDifferential transform method

${\text{EN}}_2$Error norm

ExpFDMExponential finite difference method

F-NFitzhugh–Nagumo

HBMHybrid block method

ICInitial condition

MDHBMMultiderivative hybrid block method

${\text{ME}}_{\infty }$Maximum absolute error

NFEvalNumber of function evaluations

NSFDNonstandard finite difference

OMDHBMOptimized multiderivative hybrid block method

PDQMPolynomial differential quadrature method

$\text{SCM}$Spectral collocation method

$\text{SOMDHBM}$Spectral optimized multiderivative hybrid block method

$\text{VIM}$Variational iteration method

References

[1] A. Wazwaz, Partial Differential Equations and Solitary Waves Theory, 2010.

[2] J. D. Murray, "Spatial Models and Biomedical Applications," Mathematical Biology, 2003.

[3] R. FitzHugh, "Impulses and Physiological States in Theoretical Models of Nerve Membrane," Biophysical Journal, vol. 1 no. 6, pp. 445-466, DOI: 10.1016/s0006-3495(61)86902-6, 1961.

[4] J. Nagumo, S. Arimoto, S. Yoshizawa, "An Active Pulse Transmission Line Simulating Nerve Axon," Proceedings of the IRE, vol. 50 no. 10, pp. 2061-2070, DOI: 10.1109/jrproc.1962.288235, 1962.

[5] S. Abbasbandy, "Soliton Solutions for the Fitzhugh–Nagumo Equation With the Homotopy Analysis Method," Applied Mathematical Modelling, vol. 32 no. 12, pp. 2706-2714, DOI: 10.1016/j.apm.2007.09.019, 2008.

[6] H. A. Abdusalam, "Analytic and Approximate Solutions for Nagumo Telegraph Reaction Diffusion Equation," Applied Mathematics and Computation, vol. 157 no. 2, pp. 515-522, DOI: 10.1016/j.amc.2003.08.050, 2004.

[7] D. G. Aronson, H. f. Weinberger, "Multidimensional Nonlinear Diffusion Arising in Population Genetics," Advances in Mathematics, vol. 30 no. 1, pp. 33-76, DOI: 10.1016/0001-8708(78)90130-5, 1978.

[8] P. Browne, E. Momoniat, F. M. Mahomed, "A Generalized Fitzhugh–Nagumo Equation," Nonlinear Analysis: Theory, Methods & Applications, vol. 68 no. 4, pp. 1006-1015, DOI: 10.1016/j.na.2006.12.001, 2008.

[9] S. S. A. Ramji, K. Eddy, A. Mujahid, "Uncountable Infinite Exact Solutions to the Fitzhugh-Nagumo Model," 2024.

[10] A. M. Shahid, M. Abbas, E. Kwessi, "Nonstandard Nearly Exact Analysis of the Fitzhugh–Nagumo Model," Symmetry, vol. 16 no. 5,DOI: 10.3390/sym16050585, 2024.

[11] I. Hashim, M. S. M. Noorani, B. Batiha, "A Note on the Adomian Decomposition Method for the Generalized Huxley Equation," Applied Mathematics and Computation, vol. 181 no. 2, pp. 1439-1445, DOI: 10.1016/j.amc.2006.03.011, 2006.

[12] B. Batiha, M. S. M. Noorani, I. Hashim, "Numerical Simulation of the Generalized Huxley Equation by He’s Variational Iteration Method," Applied Mathematics and Computation, vol. 186 no. 2, pp. 1322-1325, DOI: 10.1016/j.amc.2006.07.166, 2007.

[13] B. Batiha, M. S. M. Noorani, I. Hashim, "Application of Variational Iteration Method to the Generalized Burgers–Huxley Equation," Chaos, Solitons and Fractals, vol. 36 no. 3, pp. 660-663, DOI: 10.1016/j.chaos.2006.06.080, 2008.

[14] D. Uma, H. Jafari, S. Raja Balachandar, S. G. Venkatesh, S. Vaidyanathan, "An Approximate Solution for Stochastic Fitzhugh–Nagumo Partial Differential Equations Arising in Neurobiology Models," Mathematical Methods in the Applied Sciences, vol. 48 no. 3, pp. 2980-2998, DOI: 10.1002/mma.10471, 2025.

[15] J. Biazar, F. Mohammadi, "Application of Differential Transform Method to the Generalized Burgers–Huxley Equation," Applications and Applied Mathematics: International Journal, vol. 5 no. 2, 2010.

[16] G. Hariharan, K. Kannan, "Haar Wavelet Method for Solving Fitzhugh-Nagumo Equation," International Journal of Mathematical and Statistical Sciences, vol. 2 no. 2, pp. 59-63, 2010.

[17] R. A. Van Gorder, K. Vajravelu, K. Vajravelu, "A Variational Formulation of the Nagumo Reaction: Diffusion Equation and the Nagumo Telegraph Equation," Nonlinear Analysis: Real World Applications, vol. 11 no. 4, pp. 2957-2962, DOI: 10.1016/j.nonrwa.2009.10.016, 2010.

[18] G. Van, A. Robert, "Gaussian Waves in the Fitzhugh–Nagumo Equation Demonstrate One Role of the Auxiliary Function H," .

[19] M. Shih, E. Momoniat, F. M. Mahomed, "Approximate Conditional Symmetries and Approximate Solutions of the Perturbed Fitzhugh–Nagumo Equation," Journal of Mathematical Physics, vol. 46 no. 2,DOI: 10.1063/1.1839276, 2005.

[20] H. Li, Y. Guo, "New Exact Solutions to the Fitzhugh–Nagumo Equation," Applied Mathematics and Computation, vol. 180 no. 2, pp. 524-528, DOI: 10.1016/j.amc.2005.12.035, 2006.

[21] M. C. Nucci, P. A. Clarkson, "The Nonclassical Method Is More General Than the Direct Method for Symmetry Reductions. An Example of the Fitzhugh-Nagumo Equation," Physics Letters A, vol. 164 no. 1, pp. 49-56, DOI: 10.1016/0375-9601(92)90904-z, 1992.

[22] A. H. Bhrawy, "An Efficient Jacobi Pseudospectral Approximation for Nonlinear Complex Generalized Zakharov System," Applied Mathematics and Computation, vol. 247, pp. 30-46, DOI: 10.1016/j.amc.2014.08.062, 2014.

[23] A. H. Bhrawy, "A Jacobi Spectral Collocation Method for Solving Multi-Dimensional Nonlinear Fractional Sub-Diffusion Equations," Numerical Algorithms, vol. 73 no. 1, pp. 91-113, DOI: 10.1007/s11075-015-0087-2, 2016.

[24] A. H. Bhrawy, E. H. Doha, S. S. Ezz-Eldien, M. A. Abdelkawy, "A Numerical Technique Based on the Shifted Legendre Polynomials for Solving the Time-Fractional Coupled Kdv Equations," Calcolo, vol. 53,DOI: 10.1007/s10092-014-0132-x, 2016.

[25] R. E. Mickens, Advances in the Applications of Nonstandard Finite Diffference Schemes, 2005.

[26] L. I. W. Roeger, R. E. Mickens, "Exact Finite-Difference Schemes for First Order Differential Equations Having Three Distinct Fixed-Points," Journal of Difference Equations and Applications, vol. 13 no. 12, pp. 1179-1185, DOI: 10.1080/10236190701466439, 2007.

[27] L. I. W. Roeger, "Exact Finite-Difference Schemes for Two-Dimensional Linear Systems With Constant Coefficients," Journal of Computational and Applied Mathematics, vol. 219 no. 1, pp. 102-109, DOI: 10.1016/j.cam.2007.07.015, 2008.

[28] L. I. W. Roeger, "Exact Nonstandard Finite-Difference Methods for A Linear System-the Case of Centers," Journal of Difference Equations and Applications, vol. 14 no. 4, pp. 381-389, DOI: 10.1080/10236190701607669, 2008.

[29] S. Rucker, "Exact Finite Difference Scheme for an Advection-Reaction Equation," Journal of Difference Equations and Applications, vol. 9 no. 11, pp. 1007-1013, DOI: 10.1080/1023619031000146931, 2003.

[30] E. H. Doha, B. Dumitru, A. H. Bhrawy, M. A. Abdelkawy, "A Jacobi Collocation Method for Solving Nonlinear Burgers‐Type Equations," Abstract and Applied Analysis 2013, vol. 1, 2013.

[31] M. Y. Hussaini, L. S. Craig, A. Z. Thomas, "Spectral Methods for Partial Differential Equations," 1983.

[32] L. N. Trefethen, "Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations," 1996.

[33] P. J. Davis, "Interpolation and Approximation," Courier Corporation, 1975.

[34] A. H. Khater, R. S. Temsah, D. k. Callebaut, "Numerical Solutions for Some Coupled Nonlinear Evolution Equations by Using Spectral Collocation Method," Mathematical and Computer Modelling, vol. 48 no. 7-8, pp. 1237-1253, DOI: 10.1016/j.mcm.2008.02.001, 2008.

[35] R. Jiwari, "A Haar Wavelet Quasilinearization Approach for Numerical Simulation of Burgers’ Equation," Computer Physics Communications, vol. 183 no. 11, pp. 2413-2423, DOI: 10.1016/j.cpc.2012.06.009, 2012.

[36] M. Kamrani, S. M. Hosseini, "Spectral Collocation Method for Stochastic Burgers Equation Driven by Additive Noise," Mathematics and Computers in Simulation, vol. 82 no. 9, pp. 1630-1644, DOI: 10.1016/j.matcom.2012.03.007, 2012.

[37] Y. Sun, B. W. Li, "Chebyshev Collocation Spectral Method for One-Dimensional Radiative Heat Transfer in Graded Index Media," International Journal of Thermal Sciences, vol. 48 no. 4, pp. 691-698, DOI: 10.1016/j.ijthermalsci.2008.07.003, 2009.

[38] Y. Sun, B. W. Li, "Chebyshev Collocation Spectral Approach for Combined Radiation and Conduction Heat Transfer in One-Dimensional Semitransparent Medium With Graded Index," International Journal of Heat and Mass Transfer, vol. 53 no. 7-8, pp. 1491-1497, DOI: 10.1016/j.ijheatmasstransfer.2009.11.047, 2010.

[39] R. C. Mittal, R. Jiwari, "Numerical Solution of Two-Dimensional Reaction-Diffusion Brusselator System," Applied Mathematics and Computation, vol. 217 no. 12, pp. 5404-5415, DOI: 10.1016/j.amc.2010.12.010, 2011.

[40] R. C. Mittal, R. Jiwari, "A Differential Quadrature Method for Numerical Solutions of Burgers’-Type Equations," International Journal of Numerical Methods for Heat and Fluid Flow, vol. 22 no. 7, pp. 880-895, DOI: 10.1108/09615531211255761, 2012.

[41] S. S. Motsa, V. M. Magagula, P. Sibanda, "A Bivariate Chebyshev Spectral Collocation Quasilinearization Method for Nonlinear Evolution Parabolic Equations," The Scientific World Journal, vol. 2014 no. 1,DOI: 10.1155/2014/581987, 2014.

[42] A. Rezazadeh, M. Mahmoud, D. Majid, "Space–Time Spectral Collocation Method for One-Dimensional PDE Constrained Optimization," International Journal of Control, vol. 93 no. 5, pp. 1231-1241, 2020.

[43] A. Shadrin, "Error Bounds for Lagrange Interpolation," Journal of Approximation Theory, vol. 80 no. 1, pp. 25-49, DOI: 10.1006/jath.1995.1003, 1995.

[44] F. B. Khumalo, S. S. Motsa, V. M. Magagula, "A Bivariate Spectral Linear Partition Method for Solving Nonlinear Evolution Equations," Mathematical Methods in the Applied Sciences, vol. 47 no. 7, pp. 6607-6621, DOI: 10.1002/mma.9940, 2024.

[45] J. Huang, N. Nie, Y. Tang, "A Second Order Finite Difference-Spectral Method for Space Fractional Diffusion Equations," Science China Mathematics, vol. 57 no. 6, pp. 1303-1317, DOI: 10.1007/s11425-013-4716-8, 2014.

[46] M. A. Rufai, A. A. Kosti, Z. A. Anastassi, B. Carpentieri, "A New Two-Step Hybrid Block Method for the FitzHugh–Nagumo Model Equation," Mathematics (Basel), vol. 12 no. 1,DOI: 10.3390/math12010051, 2023.

[47] E. O. Adeyefa, E. O. Omole, A. Shokri, "Numerical Solution of Second-Order Nonlinear Partial Differential Equations Originating From Physical Phenomena Using Hermite Based Block Methods," Results in Physics, vol. 46,DOI: 10.1016/j.rinp.2023.106270, 2023.

[48] E. Olaoluwa Omole, E. Olusheye Adeyefa, K. Iyabo Apanpa, V. Iyadunni Ayodele, F. Emmanuel Amoyedo, H. Emadifar, "Unveiling the Power of Implicit Six-Point Block Scheme: Advancing Numerical Approximation of Two-Dimensional PDEs in Physical Systems," PLoS One, vol. 19 no. 5,DOI: 10.1371/journal.pone.0301505, 2024.

[49] A. Mehta, G. Singh, H. Ramos, "Numerical Solution of Time Dependent Nonlinear Partial Differential Equations Using a Novel Block Method Coupled With Compact Finite Difference Schemes," Computational and Applied Mathematics, vol. 42 no. 4,DOI: 10.1007/s40314-023-02345-3, 2023.

[50] M. A. Rufai, H. Ramos, "A Variable Step-Size Fourth-Derivative Hybrid Block Strategy for Integrating Third-Order IVPs, With Applications," International Journal of Computer Mathematics, vol. 99 no. 2, pp. 292-308, DOI: 10.1080/00207160.2021.1907357, 2022.

[51] U. O. Rufai, P. Sibanda, S. P. Goqo, "A One-Step Multi-Derivative Hybrid Block Method With Modified-Picard Iteration for the Solution of Second Order IVP," Engineering Letters, vol. 31 no. 4, 2023.

[52] S. J. Kayode, B. T. Olabode, O. S. Ige, "Four-Step Multi-Derivative Predictor Corrector Method for the Solution of Second Order Ordinary Differential Equations," World Scientific News, vol. 176, pp. 133-144, 2023.

[53] U. O. Rufai, P. Sibanda, S. Goqo, "An Overlapping Adaptive Step-Size Multi-Derivative Hybrid Block Method for Higher Order Initial Value Problems," Computational Methods for Differential Equations, 2025.

[54] H. Ramos, Z. Kalogiratou, T. Monovasilis, T. E. Simos, "An Optimized Two-Step Hybrid Block Method for Solving General Second Order Initial-Value Problems," Numerical Algorithms, vol. 72 no. 4, pp. 1089-1102, DOI: 10.1007/s11075-015-0081-8, 2016.

[55] H. Ramos, "An Optimized Two-Step Hybrid Block Method for Solving rst-Order Initial-Value Problems in ODEs," Differential Geometry-Dynamical Systems, vol. 19, 2017.

[56] B. S. Kashkari, S. Alqarni, "Optimization of Two-Step Block Method With Three Hybrid Points for Solving Third Order Initial Value Problems," The Journal of Nonlinear Science and Applications, vol. 12 no. 07, pp. 450-469, DOI: 10.22436/jnsa.012.07.04, 2019.

[57] S. Motsa, "Overlapping Grid-Based Optimized Single-Step Hybrid Block Method for Solving First-Order Initial Value Problems," Algorithms, vol. 15 no. 11,DOI: 10.3390/a15110427, 2022.

[58] G. Singh, A. Garg, V. Kanwar, H. Ramos, "An Efficient Optimized Adaptive Step-Size Hybrid Block Method for Integrating Differential Systems," Applied Mathematics and Computation, vol. 362,DOI: 10.1016/j.amc.2019.124567, 2019.

[59] R. Singla, G. Singh, H. Ramos, V. Kanwar, "An Efficient Optimized Adaptive Step-Size Hybrid Block Method for Integrating Directly," Journal of Computational and Applied Mathematics, vol. 420,DOI: 10.1016/j.cam.2022.114838, 2023.

[60] S. D. Yakubu, P. Sibanda, "One‐Step Family of Three Optimized Second‐Derivative Hybrid Block Methods for Solving First‐Order Stiff Problems," Journal of Applied Mathematics, vol. 2024 no. 1,DOI: 10.1155/2024/5078943, 2024.

[61] R. Jiwari, R. K. Gupta, V. Kumar, "Polynomial Differential Quadrature Method for Numerical Solutions of The Generalized Fitzhugh–Nagumo Equation with Time-Dependent Coefficients," Ain Shams Engineering Journal, vol. 5 no. 4, pp. 1343-1350, DOI: 10.1016/j.asej.2014.06.005, 2014.

[62] M. Namjoo, S. Zibaei, "Numerical Solutions of Fitzhugh–Nagumo Equation by Exact Finite-Difference and NSFD Schemes," Computational and Applied Mathematics, vol. 37 no. 2, pp. 1395-1411, DOI: 10.1007/s40314-016-0406-9, 2018.

[63] B. İnan, K. K. Ali, A. Saha, T. Ak, "Analytical and Numerical Solutions of the Fitzhugh–Nagumo Equation and Their Multistability Behavior," Numerical Methods for Partial Differential Equations, vol. 37 no. 1,DOI: 10.1002/num.22516, 2021.

[64] S. A. Manaa, F. H. Easif, A. S. Faris, "The Finite Difference Methods for FitzHugh-Nagumo Equation," Journal of Mathematics, vol. 11 no. 2, pp. 51-55, 2015.

[65] S. A. A. Ahmedai Abd Allah, P. Sibanda, S. P. Goqo, U. O. Rufai, H. Sithole Mthethwa, O. A. I. Noreldin, "A Block Hybrid Method With Equally Spaced Grid Points for Third-Order Initial Value Problems," AppliedMath, vol. 4 no. 1, pp. 320-347, DOI: 10.3390/appliedmath4010017, 2024.

[66] S. Ola Fatunla, "Block Methods for Second Order Odes," International Journal of Computer Mathematics, vol. 41 no. 1–2, pp. 55-63, DOI: 10.1080/00207169108804026, 1991.

[67] R. Singla, G. Singh, H. Ramos, V. Kanwar, "A Family of A‐Stable Optimized Hybrid Block Methods for Integrating Stiff Differential Systems," Mathematical Problems in Engineering, vol. 1, 2022.

[68] L. N. Trefethen, Spectral Methods in MATLAB, 2000.

[69] S. S. Motsa, "On the Optimal Auxiliary Linear Operator for the Spectral Homotopy Analysis Method Solution of Nonlinear Ordinary Differential Equations," Mathematical Problems in Engineering, vol. 2014 no. 1,DOI: 10.1155/2014/697845, 2014.