Numerical Approach for Solving Delay Differential Equations with Boundary Conditions

Abstract

In the present paper, a fifth-order direct multistep block method is proposed for solving the second-order Delay Differential Equations (DDEs) directly with boundary conditions using constant step size. In many life sciences applications, a delay plays an essential role in modelling natural phenomena with data simulation. Thus, an efficient numerical method is needed for the numerical treatment of time delay in the applications. The proposed direct block method computes the numerical solutions at two points concurrently at each computed step along the interval. The types of delays involved in this research are constant delay, pantograph delay, and time-dependent delay. The shooting technique is utilized to deal with the boundary conditions by applying a Newton-like method to guess the next initial values. The analysis of the proposed method based on the order, consistency, convergence, and stability of the method are discussed in detail. Four tested problems are presented to measure the efficiency of the developed direct multistep block method. The numerical simulation indicates that the proposed direct multistep block method performs better than existing methods in terms of accuracy, total function calls, and execution times.

1. Introduction

Various problems in sciences and engineering may be interpreted as mathematical models in terms of Delay Differential Equations (DDEs). DDEs exist in many real-life problems, such as in dengue fever diseases where time delay occurs between the mosquito bite until the body gets the infection. There is also an example in the biological process where delay exists during cell division and driver response’s time in transport delay. Bocharov et al. [1] employed DDEs and assimilation with the models of real data in studying the immune system and dynamic of infections. Thus, a time delay is naturally present in most problems, which makes it impossible to be ignored [2].

The general second order DDEs are given as:

$$y^{″}=f\left(x,y\left(x\right),y^{\prime }\left(x\right),y(x-\tau ),y^{\prime }(x-\tau )\right),x\in [a,b]$$

(1)

$$y\left(x\right)=\varphi \left(x\right),x\in [a-\tau ,a],\tau \in {\mathcal{\Re }}^{+}$$

(2)

subject to boundary conditions:

$$y\left(a\right)=\alpha ,y\left(b\right)=\beta ,$$

(3)

where $\alpha =\varphi \left(a\right)$, $\tau$ is a positive constant delay term for constant delay type, or a non-constant delay term, which is a function of x, $\tau =\tau \left(x\right)$ for pantograph and time-dependent delay type while $\varphi \left(x\right)$ is the smooth initial function given in the problems. For constant delay type, $\tau =mh$, where $m=1,2,3,\dots$ and h is the step size. The general form of pantograph delay type is shown as follows:

$$y^{″}=f\left(x,y\left(x\right),y^{\prime }\left(x\right),y\left(qx\right),y^{\prime }\left(qx\right)\right),x\in [a,b]$$

(4)

where $0<q<1$.

The presence of the delay term makes the solutions of DDEs fall in the previous time giving apparent distinction from Ordinary Differential Equations (ODEs) since the solutions of ODEs fall in the current time.

The study on solving the second-order DDEs of constant delay type numerically with boundary conditions started in the 1970s, where Nevers and Schmitt [3] solved this problem by using Euler’s method and shooting technique. Other numerical methods such as finite differences method, cubic splines method, midpoint difference method, collocation method and Richardson extrapolation method have been used to solve second-order DDEs in [4,5,6,7,8]. However, only a few studies have numerically solved the second-order DDEs of non-constant delay type with boundary conditions. For example, Bakke and Jackiewicz [8] acquired the convergent method of order three when applying Richardson extrapolation while Cahlon and Nachman [9] implemented Atkinson’s product integration technique to solve time-dependent delay problems. Apart from that, Bica et al. [10] implemented the combination of the Picard sequence and the trapezoidal rule to solve four pantograph delay problems, approximated by using the piecewise Birkhoff cubic interpolation. Besides that, Wazwaz et al. [11] computed four nonlinear and linear pantograph delay problems by utilizing the analytical methods, Adomian decomposition, and the variational iteration method.

The direct approach for solving second-order DDEs with initial conditions using a direct multistep method has been proposed by Seong et al. [12,13,14]. The proposed direct methods in [12,13] were implemented using the predictor–corrector technique to calculate the one-point approximate solution for each iteration. Then, Seong and Majid [14] presented a direct two-point block method to compute two approximate solutions simultaneously at two different points. The approach used in [12,13,14] is in a fully implicit manner.

The implementation of the shooting technique has been discussed in [15,16,17,18] to solve the second order Boundary Value Problems (BVPs) using direct multistep methods. Majid et al. [15] and Phang et al. [16] have implemented the Newton’s method to obtain the missing initial values for solving BVPs of Dirichlet and Neumann boundary conditions, respectively. Phang et al. [17] has adapted the multiple shooting techniques and employed the three-step iterative method to produce the missing guessing values when solving BVPs problems with both Dirichlet and Neumann boundary conditions. Meanwhile, Mohd Nasir et al. [18] had solved the boundary value problems with a Robin boundary condition using Steffensen’s method to guess the missing initial values in shooting technique.

This paper aims to propose a direct multistep block method in the form of predictor corrector block method for solving the second-order DDEs with boundary conditions directly. The shooting technique has been adapted to obtain the missing initial value, while the Newton-like method proposed in [19] will be used to determine the next initial values. The numerical computations have been performed using the C language.

2. Derivation of Method

The interval [$a,b$] is partitioned into a series of a block, as portrayed in Figure 1, where h is the step size.

This study will consider a two-point block approach to approximate two solutions for $y_{i+1}$ and $y_{i+2}$ at the points $x_{i+1}$ and $x_{i+2}$, respectively, at each step. The backward points taken for every block of the method of order five are selected in a diagonally implicit way. The backward points are {$x_i,x_{i-1},x_{i-2},x_{i-3}$} and {$x_{i+1},x_i,x_{i-1},x_{i-2}$} for the points $y_{i+1}$ and $y_{i+2}$, respectively.

Next, the derivation formulae for the first point, $y_{i+1}$ and second point, $y_{i+2}$ are obtained by taking the integration once and twice on both sides for the second order ODEs defined as follows:

$$y^{″}=f\left(x,y,y^{\prime }\right).$$

First point, $y_{i+1}$:

Integrate once:

$${\int }_{x_i}^{x_{i+1}}{y}^{″}dx={\int }_{x_i}^{x_{i+1}}f\left(x,y,y^{\prime }\right)dx$$

$$y_{i+1}^{\prime }=y_i^{\prime }+{\int }_{x_i}^{x_{i+1}}f\left(x,y,y^{\prime }\right)dx$$

(5)

Integrate twice:

$${\int }_{x_i}^{x_{i+1}}{\int }_{x_i}^x{y}^{″}dxdx={\int }_{x_i}^{x_{i+1}}{\int }_{x_i}^{x}f\left(x,y,y^{\prime }\right)dxdx$$

$$y_{i+1}=y_i+h{y}_i^{\prime }+{\int }_{x_i}^{x_{i+1}}(x_{i+1}-x)f\left(x,y,y^{\prime }\right)dx,$$

(6)

where $h=x_{i+1}-x_i.$

Second point, $y_{i+2}$:

Integrate once:

$${\int }_{x_i}^{x_{i+2}}{y}^{″}dx={\int }_{x_i}^{x_{i+2}}f\left(x,y,y^{\prime }\right)dx$$

$$y_{i+2}^{\prime }=y_i^{\prime }+{\int }_{x_i}^{x_{i+2}}f\left(x,y,y^{\prime }\right)dx.$$

(7)

Integrate twice:

$${\int }_{x_i}^{x_{i+2}}{\int }_{x_i}^x{y}^{″}dxdx={\int }_{x_i}^{x_{i+2}}{\int }_{x_i}^{x}f\left(x,y,y^{\prime }\right)dxdx$$

$$y_{i+2}=y_i+2h{y}_i^{\prime }+{\int }_{x_i}^{x_{i+2}}(x_{i+2}-x)f\left(x,y,y^{\prime }\right)dx.$$

(8)

Then, the function $f\left(x,y,y^{\prime }\right)$ in Equations (5)–(8) will be approximated by using Lagrange interpolation [18]:

$$P_q\left(x\right)=\sum _{j=0}^q\left(\prod _{n=0,n\ne j}^q\frac{x-x_{i+p-n}}{x_{i+p-j}-x_{i+p-n}}\right)f_{i+p-j},p=0,1,2,$$

where $P_q\left(x\right)$ is the Lagrange polynomial of degree q. The proposed direct multistep block method will consists of predictor and corrector block methods. For the predictor formula, taking $p=0$ for both first and second point while the degree of Lagrange polynomial is $q=3$. The points involved to interpolate the predictor for both $y_{i+1}$ and $y_{i+2}$ is as shown below.

$$\{(x_i,f_i),(x_{i-1},f_{i-1}),(x_{i-2},f_{i-2}),(x_{i-3},f_{i-3})\}.$$

For the corrector formula, taking $p=1$ for the first point and $p=2$ for the second point while the degree of Lagrange polynomial is $q=4$. The points involved to interpolate the corrector for the first point, $y_{i+1}$ is as shown below

$$\{(x_{i+1},f_{i+1}),(x_i,f_i),(x_{i-1},f_{i-1}),(x_{i-2},f_{i-2}),(x_{i-3},f_{i-3})\}$$

and for the second point, $y_{i+2}$ is as follows

$$\{(x_{i+2},f_{i+2}),(x_{i+1},f_{i+1}),(x_i,f_i),(x_{i-1},f_{i-1}),(x_{i-2},f_{i-2})\}.$$

By taking $x=x_{i+2}+sh$ and $dx=hds$, then substitute into the integral part. By evaluating and simplifying the integration limit from −2 to 1 and −2 to 0 of the first and second point, respectively using MAPLE. Hence, the corrector formula of the first and second point of the direct multistep block method (2PBM5) is obtained as follows:

First point, $y_{i+1}$:

$$\begin{array}{cc}\hfill y_{i+1}^{\prime }& =y_i^{\prime }+\frac{h}{720}\left(251f_{i+1}+646f_i-264f_{i-1}+106f_{i-2}-19f_{i-3}\right),\hfill \\ \hfill y_{i+1}& =y_i+h{y}_i^{\prime }+\frac{h^2}{1440}\left(135f_{i+1}+752f_i-246f_{i-1}+96f_{i-2}-17f_{i-3}\right),\hfill \end{array}$$

(9)

Second point, $y_{i+2}$:

$$\begin{array}{cc}\hfill y_{i+2}^{\prime }& =y_i^{\prime }+\frac{h}{90}\left(29f_{i+2}+124f_{i+1}+24f_i+4f_{i-1}-f_{i-2}\right),\hfill \\ \hfill y_{i+2}& =y_i+2h{y}_i^{\prime }+\frac{h^2}{90}\left(5f_{i+2}+104f_{i+1}+78f_i-8f_{i-1}+f_{i-2}\right).\hfill \end{array}$$

(10)

The proposed method will be implemented with the predictor–corrector technique, where the corrector method possesses one order higher than the predictor method. The development of the predictor formula can be obtained using the same approach as the corrector formula. Thus, the predictor formula for the first and second point is as shown below:

First point, $y_{i+1}$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{i+1}^{\prime }& =y_i^{\prime }+\frac{h}{24}\left(55f_i-59f_{i-1}+37f_{i-2}-9f_{i-3}\right),\hfill \\ \hfill y_{i+1}& =y_i+h{y}_i^{\prime }+\frac{h^2}{360}\left(323f_i-264f_{i-1}+159f_{i-2}-38f_{i-3}\right),\hfill \end{array}\end{array}$$

Second point, $y_{i+2}$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{i+2}^{\prime }& =y_i^{\prime }+\frac{h}{3}\left(27f_i-44f_{i-1}+31f_{i-2}-8f_{i-3}\right),\hfill \\ \hfill y_{i+2}& =y_i+2h{y}_i^{\prime }+\frac{h^2}{45}\left(272f_i-366f_{i-1}+246f_{i-2}-62f_{i-3}\right).\hfill \end{array}\end{array}$$

3. Analysis of Method

3.1. Order of Method

The derived 2PBM5 method could be interpreted by employing the matrix difference equation, as shown below

$$\alpha Y_M=h\beta Y_M^{\prime }+h^2\gamma F_M,$$

where

$$Y_M=\left[\begin{array}{c}{y}_{i-2}\\ y_{i-1}\\ y_i\\ y_{i+1}\\ y_{i+2}\end{array}\right],Y_M^{\prime }=\left[\begin{array}{c}{y}_{i-2}^{\prime }\\ y_{i-1}^{\prime }\\ y_i^{\prime }\\ y_{i+1}^{\prime }\\ y_{i+2}^{\prime }\end{array}\right],F_M=\left[\begin{array}{c}{f}_{i-3}\\ f_{i-2}\\ f_{i-1}\\ f_i\\ f_{i+1}\\ f_{i+2}\end{array}\right].$$

The linear difference operator as defined in Fatunla [20] is given by:

$$L[y\left(x\right);h]=\sum _{j=0}^k[{\alpha }_{j}y(x+jh)-h{\beta }_j{y}^{\prime }(x+jh)-h^2{\gamma }_j{y}^{″}(x+jh)],$$

(11)

where $y\left(x\right)$ is a smooth function on $[a,b]$. The dependent variable and its derivatives are substituted by using the Taylor expansion. By gathering the terms of Equation (11) eventually, we obtain the following

$$L[y\left(x\right);h]=C_{0}y\left(x\right)+C_{1}h{y}^{\left(1\right)}\left(x\right)+\dots +C_r{h}^r{y}^{\left(r\right)}\left(x\right)+\dots ,$$

where unknowns $C_r$ are in terms of constants, given as

$$\begin{array}{cc}\hfill C_0& ={\alpha }_0+{\alpha }_1+{\alpha }_2+\dots +{\alpha }_s,\hfill \\ \hfill C_1& ={\alpha }_1+2{\alpha }_2+3{\alpha }_3+\dots +s{\alpha }_s-({\beta }_0+{\beta }_1+{\beta }_2+\dots +{\beta }_s),\hfill \\ \hfill C_2& ={\displaystyle \frac{1}{2!}}({\alpha }_1+2^2{\alpha }_2+3^2{\alpha }_3+\dots +s{\alpha }_s)-({\beta }_1+2{\beta }_2+3{\beta }_3+\dots +s{\beta }_s)-({\gamma }_0+{\gamma }_1+{\gamma }_2+\dots +s^{q-2}{\gamma }_s),\hfill \\ \hfill ⋮\\ \hfill C_r& =\frac{1}{r!}({\alpha }_1+2^r{\alpha }_2+3^r{\alpha }_3+\dots +s^r{\alpha }_s)-{\displaystyle \frac{1}{(r-1)!}}({\beta }_1+2^{r-1}{\beta }_2+3^{r-1}{\beta }_3+\dots +s^{r-1}{\beta }_s)\hfill \\ \hfill & -\frac{1}{(r-2)!}({\gamma }_1+2^{r-2}{\gamma }_2+3^{r-2}{\gamma }_3+\dots +s^{r-2}{\gamma }_s),\hfill \\ \hfill & r=3,4,5,\dots \hfill \end{array}$$

(12)

Definition 1.

(Fatunla [20]) The order r of the linear difference operator L[y(x);h] is a unique integer r such that $C_{r+2}\ne 0$.

The non-zero column vector is the error constant of the method.

Hence, the 2PBM5 in Equations (9) and (10) can be written in the following matrix difference structure

$$\begin{array}{cc}\hfill \left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& -1& 1& 0\\ 0& 0& 0& 0& 0\\ 0& 0& -1& 0& 1\end{array}\right]\left[\begin{array}{c}{y}_{i-2}\\ y_{i-1}\\ y_i\\ y_{i+1}\\ y_{i+2}\end{array}\right]& =h\left[\begin{array}{ccccc}0& 0& 1& -1& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 1& 0& -1\\ 0& 0& 2& 0& 0\end{array}\right]\left[\begin{array}{c}{y}_{i-2}^{\prime }\\ y_{i-1}^{\prime }\\ y_i^{\prime }\\ y_{i+1}^{\prime }\\ y_{i+2}^{\prime }\end{array}\right]\hfill \\ \hfill & +h^2\left[\begin{array}{cccccc}-{\displaystyle \frac{19}{720}}& {\displaystyle \frac{53}{360}}& -{\displaystyle \frac{11}{30}}& {\displaystyle \frac{323}{360}}& {\displaystyle \frac{251}{720}}& 0\\ & & \\ -{\displaystyle \frac{17}{1440}}& {\displaystyle \frac{1}{15}}& -{\displaystyle \frac{41}{240}}& {\displaystyle \frac{47}{90}}& {\displaystyle \frac{3}{32}}& 0\\ & & \\ 0& -{\displaystyle \frac{1}{90}}& {\displaystyle \frac{2}{45}}& {\displaystyle \frac{4}{15}}& {\displaystyle \frac{62}{45}}& {\displaystyle \frac{29}{90}}\\ & & \\ 0& {\displaystyle \frac{1}{90}}& -{\displaystyle \frac{4}{45}}& {\displaystyle \frac{13}{15}}& {\displaystyle \frac{52}{45}}& {\displaystyle \frac{1}{18}}\end{array}\right]\left[\begin{array}{c}{f}_{i-3}\\ f_{i-2}\\ f_{i-1}\\ f_i\\ f_{i+1}\\ f_{i+2}\end{array}\right].\hfill \end{array}$$

(13)

By applying the associated set of coefficients from Equation (12), we eventually obtain $C_0=C_1=C_2=C_3=\cdots =C_6={\left[\begin{array}{cccc}0& 0& 0& 0\end{array}\right]}^T$ while $C_7$ is as follows

$$\begin{array}{cc}\hfill C_7& ={\left[\begin{array}{cccc}-\frac{3}{160}& -\frac{41}{5040}& -\frac{1}{90}& \frac{1}{315}\end{array}\right]}^T\ne {\left[\begin{array}{cccc}0& 0& 0& 0\end{array}\right]}^T.\hfill \end{array}$$

(14)

By referring to Definition 1, the method achieves order r if $C_0=C_1=C_2=\dots =C_{r+1}=0$ and $C_{r+2}\ne 0$ with error constant of $C_{r+2}$. When $C_7\ne 0,$ $r=5$ is obtained, which indicates that the method is of order five.

3.2. Zero Stability

Zero stability only investigates the stability of the method in the limit of the step size, $h\to 0$, while the stability theory examines the stability of the method when h is a fixed non-zero value.

The linear multistep method is defined to be zero-stable when the first characteristic polynomial is described as $\rho \left(R\right)=det\left[{\sum }_{i=0}^p{A}_{\left(i\right)}{R}^{p-i}\right]=0$, have roots $R_j$ that satisfies $|R_j|\le 1$.

The corrector formula for 2PBM5 in Equations (9) and (10) can be illustrated in the matrix form as follows:

$$\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{y}_{i+1}^{\prime }\\ y_{i+1}\\ y_{i+2}^{\prime }\\ y_{i+2}\end{array}\right]=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{y}_{i-1}^{\prime }\\ y_{i-1}\\ y_i^{\prime }\\ y_i\end{array}\right]+h\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 2& 0\end{array}\right]\left[\begin{array}{c}{y}_{i-1}^{\prime }\\ y_{i-1}\\ y_i^{\prime }\\ y_i\end{array}\right]$$

$$+h\left[\begin{array}{cccc}0& 0& -\frac{19}{720}& \frac{53}{360}\\ 0& 0& 0& 0\\ 0& 0& 0& -\frac{1}{90}\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{f}_{i-5}\\ f_{i-4}\\ f_{i-3}\\ f_{i-2}\end{array}\right]+h\left[\begin{array}{cccc}-\frac{11}{30}& \frac{323}{360}& \frac{251}{720}& 0\\ 0& 0& 0& 0\\ \frac{2}{45}& \frac{4}{15}& \frac{62}{45}& \frac{29}{90}\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{f}_{i-1}\\ f_i\\ f_{i+1}\\ f_{i+2}\end{array}\right]$$

$$+h^2\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& -\frac{17}{1440}& \frac{1}{15}\\ 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{90}\end{array}\right]\left[\begin{array}{c}{f}_{i-5}\\ f_{i-4}\\ f_{i-3}\\ f_{i-2}\end{array}\right]+h^2\left[\begin{array}{cccc}0& 0& 0& 0\\ -\frac{41}{240}& \frac{47}{90}& \frac{3}{32}& 0\\ 0& 0& 0& 0\\ -\frac{4}{45}& \frac{13}{15}& \frac{52}{45}& \frac{1}{18}\end{array}\right]\left[\begin{array}{c}{f}_{i-1}\\ f_i\\ f_{i+1}\\ f_{i+2}\end{array}\right].$$

Thus, the first characteristic polynomial of 2PBM5 is

$$\rho \left(R\right)=det\left[\sum _{i=0}^p{A}_{\left(i\right)}{R}^{p-i}\right]=det\left[A_0{R}^1-A_1\right]=0,$$

where $A_0=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$ and $A_1=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$,

$$\rho \left(R\right)=det\left[\begin{array}{cccc}R& 0& -1& 0\\ 0& R& 0& -1\\ 0& 0& R-1& 0\\ 0& 0& 0& R-1\end{array}\right]=0,$$

$$R^2{(R-1)}^2=0,$$

$$R=0,0,1,1.$$

Since $|R_j|\le 1$, the proposed method achieved zero stable.

3.3. Consistency and Convergence

Based on the definition of consistency for the block method in [20,21], we can deduce a definition as the following

Definition 2.

The block method in Equation (13) is consistent as it’s order $r\ge 1$.

Therefore, since the proposed method, 2PBM5 has order five, $r=5\ge 1$ thus 2PBM5 is consistent.

Definition 3.

(Lambert [21]) The necessary and suffcient conditions for linear multistep method to be convergent are that it must be consistent and zero-stable.

Since it has been proved previously that the proposed method is consistent and zero-stable, thus the method is convergent.

3.4. Stability Region

The following linear test equation for the second order DDEs of constant delay type

$$f=ay\left(x\right)+by(x-\tau )$$

(15)

is used to obtain the stability polynomial for 2PBM5 method where $\tau =mh$.

The matrix form of 2PBM5 formula in Equations (9) and (10) is given by

$$\begin{array}{cc}\hfill A_0{Y}_{M+1}& =A_1{Y}_M+h{B}_1{Y}_M+h\sum _{j=0}^1{C}_{j+1}{F}_{M+j}+h^2\sum _{j=0}^1{D}_{j+1}{F}_{M+j},\hfill \\ \hfill A_0{Y}_{M+1}& =A_1{Y}_M+h{B}_1{Y}_M+h{C}_1{F}_M+h{C}_2{F}_{M+1}+h^2{D}_1{F}_M+h^2{D}_2{F}_{M+1},\hfill \end{array}$$

(16)

where

$$Y_{M+1}=\left[\begin{array}{c}{y}_{i+1}^{\prime }\\ y_{i+1}\\ y_{i+2}^{\prime }\\ y_{i+2}\end{array}\right],Y_M=\left[\begin{array}{c}{y}_{i-1}^{\prime }\\ y_{i-1}\\ y_i^{\prime }\\ y_i\end{array}\right],F_M=\left[\begin{array}{c}{f}_{i-5}\\ f_{i-4}\\ f_{i-3}\\ f_{i-2}\end{array}\right],F_{M+1}=\left[\begin{array}{c}{f}_{i-1}\\ f_i\\ f_{i+1}\\ f_{i+2}\end{array}\right],$$

$$A_0=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],A_1=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],B_1=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 2& 0\end{array}\right],C_1=\left[\begin{array}{cccc}0& 0& -\frac{19}{720}& \frac{53}{360}\\ 0& 0& 0& 0\\ 0& 0& 0& -\frac{1}{90}\\ 0& 0& 0& 0\end{array}\right],$$

$$C_2=\left[\begin{array}{cccc}-\frac{11}{30}& \frac{323}{360}& \frac{251}{720}& 0\\ 0& 0& 0& 0\\ \frac{2}{45}& \frac{4}{15}& \frac{62}{45}& \frac{29}{90}\\ 0& 0& 0& 0\end{array}\right],D_1=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& -\frac{17}{1440}& \frac{1}{15}\\ 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{90}\end{array}\right],D_2=\left[\begin{array}{cccc}0& 0& 0& 0\\ -\frac{41}{240}& \frac{47}{90}& \frac{3}{32}& 0\\ 0& 0& 0& 0\\ -\frac{4}{45}& \frac{13}{15}& \frac{52}{45}& \frac{1}{18}\end{array}\right].$$

Substituting the test Equation (15) into (16) gives [14]

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A_0{Y}_{M+1}=& A_1{Y}_M+h{B}_1{Y}_M+h{C}_1(a{Y}_M+b{Y}_{M-m})+h{C}_2(a{Y}_{M+1}+b{Y}_{M+1-m})+\hfill \\ \hfill & h^2{D}_1(a{Y}_M+b{Y}_{M-m})+h^2{D}_2(a{Y}_{M+1}+b{Y}_{M+1-m}),\hfill \end{array}\end{array}$$

where $m=\frac{\tau }{h}$, $m=1,2,3,\dots$. Rearranging the equation yields

$$(A_0-ha{C}_2-D_2)Y_{M+1}-(A_1+h{B}_1+ha{C}_1+D_1+hb{C}_2)Y_M-\left(hb{C}_1\right)Y_{M-m}=0.$$

(17)

Let $H_1=ha$, $H_2=hb$ and $m=1.$ The variable m is taken to be 1 because based on [22], the stability region is the largest when $m=1$ compared to other values of m.

Thus, Equation (17) becomes

$$(A_0-H_1{C}_2-D_2)t^2-(A_1+H_1{B}_1+H_1{C}_1+D_1+H_2{C}_2)t-H_2{C}_1=0.$$

(18)

The determinant of Equation (18) is defined as the following

$$\pi (H_1,H_2;t)=det\left((A_0-H_1{C}_2-D_2)t^2-(A_1+H_1{B}_1+H_1{C}_1+D_1+H_2{C}_2)t-H_2{C}_1\right)$$

The stability polynomial obtained after computing the determinant is as follows

$$\begin{array}{cc}\hfill & t^8\left(\frac{731}{1620}-\frac{12577951}{23328000}{H}_1-\frac{75537161}{139968000}{H}_1^2\right)+\hfill \\ \hfill & t^7\left(-\frac{3764}{2025}-\frac{75537161}{69984000}{H}_1{H}_2+\frac{10093477}{11664000}{H}_1-\frac{14354737}{69984000}{H}_1^2-\frac{325331}{777600}{H}_1^3-\frac{12577951}{23328000}{H}_2\right)+\hfill \\ \hfill & t^6(\frac{11401}{8100}-\frac{325331}{388800}{H}_1^2{H}_2+\frac{15003983}{34992000}{H}_1{H}_2-\frac{22789}{155520}{H}_1^3+\frac{10093477}{11664000}{H}_2-\frac{75537161}{139968000}{H}_2^2-\hfill \\ \hfill & \frac{39418517}{139968000}{H}_1^2+\frac{13864277}{23328000}{H}_1)+\hfill \\ \hfill & t^5\left(-\frac{22789}{77760}{H}_1^2{H}_2-\frac{325331}{777600}{H}_2^2{H}_1-\frac{19892117}{69984000}{H}_1{H}_2+\frac{44362703}{69984000}{H}_2^2+\frac{13864277}{23328000}{H}_2\right)+\hfill \\ \hfill & t^4\left(-\frac{22789}{155520}{H}_2^2{H}_1-\frac{365717}{139968000}{H}_2^2\right)=0.\hfill \end{array}$$

(19)

Next, the stability region is illustrated in the $(H_1-H_2)$ plane by applying the boundary locus technique and substituting $t=1,-1$ and $t=$ cos $\theta +i$ sin $\theta$, $0\le \theta \le 2\pi$ in the stability polynomial (19). The real and imaginary parts of t = $cos$ $\theta +i$ sin $\theta$ are separated and solved simultaneously to obtain the points in the region. Since the set of all roots in the stability polynomial (19) satisfies $\left|t\right|\le 1$ and lies within the region’s boundary, the stability region is absolutely stable. The shaded area in Figure 2 is the stability region for 2PBM5.

4. Implementation of Method

The implementation of the developed 2PBM5 method consists of two parts i.e., delay differential equations and boundary value problems.

4.1. Delay Differential Equation

The implementation of second order DDEs using 2PBM5 started by using Direct Euler’s method as the predictor and Direct Modified Euler’s method as the corrector to compute the three starting values, $y_1$,$y_2$ and $y_3$. After the starting values are adequate, 2PBM5 will be implemented until the end of interval, N. The location of the delay arguments for constant delay case, $\left(x-\tau \right)$, pantograph delay case, $\left(qx\right)$ or time-dependent delay case, $\left(x-\tau \left(x\right)\right)$ are stored. The implementation for constant and non-constant delay cases are as follows:

4.1.1. Constant Delay Case

If $\left(x-\tau \right)\le a$, use the initial function, $\varphi \left(x\right)$ to calculate $y\left(x-\tau \right)$ and $y^{\prime }\left(x-\tau \right)$. Otherwise, the delay solutions depend on the location of $\left(x-\tau \right)$. The location may allow us to recall the previously calculated $y\left(x-\tau \right)$ and $y^{\prime }\left(x-\tau \right)$ which had been stored earlier since the implementation in fixed step size i.e., $\tau =mh$ for $m=1,2,3,\dots .$

4.1.2. Non-Constant Delay Case (Pantograph Delay)

In the code, we implemented the simple linear interpolation to calculate the first initial delay solution, then the quadratic interpolation will be used to calculate the second delay solution for the initial second point. Later, as the computation continued, more points were stored and this allowed the code to use those points to interpolate the delay solutions using Lagrange interpolations at higher order. For each delay solution, we need to locate the delay argument occurrence, then we select the closest points to interpolate it. After the delay solution been calculated, it needs to be updated.

4.1.3. Non-Constant Delay Case (Time-Dependent Delay)

The computation strategy for the delay solutions is similar to the constant delay case and pantograph delay case, where for $\left(x-\tau \left(x\right)\right)\le a$, the delay solution is calculated using the initial function, $\varphi \left(x\right)$ to approximate $y\left(x-\tau \left(x\right)\right)$ and $y^{\prime }\left(x-\tau \left(x\right)\right)$. For $\left(x-\tau \left(x\right)\right)>a$, the delay solution depends on the location of $\left(x-\tau \left(x\right)\right)$ and we used Lagrange interpolation to interpolate the delay solution. The choices of order for Lagrange interpolations depend on the number of points that have been stored and the code will select the closest point to interpolate.

A procedure to solve the second order DDEs by using 2PBM5 can be observed in Algorithm 1 below:

Algorithm 1
Step 1:	Set the starting value a, ending value b, step size h, given boundary values and given initial function, $\varphi \left(x\right)$.
Step 2:	For $i=0,1,2,3$. Calculate $x_{i+1}=a+(i+1)h$. Compute function $f_i$ and delay argument of the equation, $(x-\tau )$, $\left(qx\right)$ or $\left(x-\tau \left(x\right)\right)$.
Step 3:	Calculate the starting values, $y_1,y_2$ and $y_3$ using Direct Euler’s method as the predictor and Direct Modified Euler’s method as the corrector.
Step 4:	For $i=4,\dots ,N$, calculate $x_{i+1}=x_i+(i+1)h$. Compute function $f_i$ and delay argument of the equation, $(x-\tau )$, $\left(qx\right)$ or $\left(x-\tau \left(x\right)\right)$.
Step 5:	Calculate $y_{i+1}^{\prime }$, $y_{i+1}$, $y_{i+2}^{\prime }$ and $y_{i+2}$ using predictor formula:
	$y_{i+1}^{\prime }=y_i^{\prime }+\frac{h}{24}\left(55f_i-59f_{i-1}+37f_{i-2}-9f_{i-3}\right),$
	$y_{i+1}=y_i+h{y}_i^{\prime }+\frac{h^2}{360}\left(323f_i-264f_{i-1}+159f_{i-2}-38f_{i-3}\right),$
	$y_{i+2}^{\prime }=y_i^{\prime }+\frac{h}{3}\left(27f_i-44f_{i-1}+31f_{i-2}-8f_{i-3}\right),$
	$y_{i+2}=y_i+2h{y}_i^{\prime }+\frac{h^2}{45}\left(272f_i-366f_{i-1}+246f_{i-2}-62f_{i-3}\right).$
	Evaluate $f_{i+1}$ and $f_{i+2}$. Calculate $y_{i+1}^{\prime }$, $y_{i+1}$, $y_{i+2}^{\prime }$ and $y_{i+2}$ using corrector formula:
	$y_{i+1}^{\prime }=y_i^{\prime }+\frac{h}{720}\left(251f_{i+1}+646f_i-264f_{i-1}+106f_{i-2}-19f_{i-3}\right),$
	$y_{i+1}=y_i+h{y}_i^{\prime }+\frac{h^2}{1440}\left(135f_{i+1}+752f_i-246f_{i-1}+96f_{i-2}-17f_{i-3}\right),$
	$y_{i+2}^{\prime }=y_i^{\prime }+\frac{h}{90}\left(29f_{i+2}+124f_{i+1}+24f_i+4f_{i-1}-f_{i-2}\right),$
	$y_{i+2}=y_i+2h{y}_i^{\prime }+\frac{h^2}{90}\left(5f_{i+2}+104f_{i+1}+78f_i-8f_{i-1}+f_{i-2}\right).$
Step 6:	Stop.

4.2. Boundary Value Problems

The shooting technique was utilized to manage the Boundary Value Problems (BVPs) engaged in this study. The BVPs were transformed to Initial Value Problems (IVPs) to generate the second initial value $y^{\prime }\left(a\right)$ with the implementation of Newton-like method in [19] since the initial value given by $y\left(a\right)$ is not sufficient. The idea of shooting technique is to guess the $y^{\prime }\left(a\right)$ as accurate as possible with less number of guessing times.

Suppose that second-order ODEs are defined such that the dependent variable, y depends on both x and parameter t as follows:

$$y^{″}(x,t)=f\left(x,y(x,t),y^{\prime }(x,t)\right),$$

(20)

with initial conditions:

$$y(a,t)=\alpha ,y^{\prime }(a,t)=t_1.$$

The parameters $t=t_k$ is chosen such that

$$\underset{k\to \infty }{lim}y(b,t_k)-\beta =0.$$

The first initial guessing, $t_1$ is obtained by using:

$$t_1={\displaystyle \frac{\beta -\alpha }{b-a}}.$$

Then, partial differentiation of Equation (20) with respect to t is computed by considering

$$z(x,t)={\displaystyle \frac{\delta y}{\delta t}}(x,t).$$

Thus, we obtain

$$z^{″}(x,t)={\displaystyle \frac{\delta f}{\delta y}}(x,y,y^{\prime })z(x,t)+{\displaystyle \frac{\delta f}{\delta y^{\prime }}}(x,y,y^{\prime })z^{\prime }(x,t),$$

(21)

with initial conditions:

$$z(a,t)=0,z^{\prime }(a,t)=1.$$

Consequently, both IVPs, Equations (20) and (21) are computed simultaneously to acquire two solutions $y(x,t)$ and $z(x,t)$ to be substituted into the Newton’s-like method in obtaining the next guessing initial value, $t_k$. The Newton’s-like method is defined as:

$$t_k=w_{k-1}-{\displaystyle \frac{y(b,w_{k-1})-\beta }{z(b,w_{k-1})}},$$

where

$$w_{k-1}=t_{k-1}-{\displaystyle \frac{y(b,t_{k-1})-\beta }{z(b,t_{k-1})}}.$$

The procedure is stopped as the absolute error, $\left|y\right(b,t_{k-1})-\beta |\le TOL$ where $TOL$ is the tolerance selected to be ${10}^{-5}$.

5. Numerical Simulation

There are four tested problems in this numerical simulation, solved by using the proposed method, 2PBM5. The two-point block fully implicit method of order 5 (2PBFM5) in [14] will be implemented using the same strategy discussed in Section 4. The results obtained are also compared with the cubic splines method [5] for constant delay problem, fixed point technique [10] for pantograph delay problems, and iterative method [9] for time-dependent delay problem.

Problem 1.

DDE of constant delay (constant delay term, $\tau =0.1$):

$y^{″}\left(x\right)=-y(x-0.1)+10$ sin$(10x-1)-1000$ sin$\left(10x\right),0\le x\le 1$

$y\left(x\right)=10$ sin$\left(10x\right),x\le 0$

$y\left(1\right)=10$ sin$\left(10\right),$

Source: Sakai [5].

Problem 2.

DDE of pantograph delay (non-constant delay term, $\tau \left(x\right)=\frac{x}{2}$):

$y^{″}\left(x\right)=-2e^{-x}+{\displaystyle \frac{y\left(x\right)}{2}}+e^{-\frac{x}{2}}y\left({\displaystyle \frac{x}{2}}\right),0\le x\le 1$

$y\left(0\right)=0,$

$y\left(1\right)=e^{-1},$

Exact solution: $y\left(x\right)=x{e}^{-x}.$

Source: Bica et al. [10].

Problem 3.

DDE of pantograph delay (non-constant delay term, $\tau \left(x\right)=\frac{x}{2}$):

$y^{″}\left(x\right)=1+2\left(1+{\displaystyle \frac{x^2}{8}}\right)$ cos$\left({\displaystyle \frac{x}{2}}\right)-2$ cos$\left({\displaystyle \frac{x}{2}}\right)y\left({\displaystyle \frac{x}{2}}\right),0\le x\le {\displaystyle \frac{\pi }{4}}$

$y\left(0\right)=1,$

$y\left({\displaystyle \frac{\pi }{4}}\right)=1+{\displaystyle \frac{\sqrt{2}}{2}}+{\displaystyle \frac{{\pi }^2}{32}},$

Exact solution: $y\left(x\right)=1+{\displaystyle \frac{x^2}{2}}+$ sin$\left(x\right).$

Source: Bica et al. [10].

Problem 4.

DDE of time-dependent delay (non-constant delay term, $\tau \left(x\right)=\frac{x}{2}-0.1$):

$y^{″}\left(x\right)=-2y^{\prime }+2y\left({\displaystyle \frac{x}{2}}-0.1\right)-2e^{-2(\frac{x}{2}-0.1)}-4\left({\displaystyle \frac{x}{2}}-0.1\right)+4,x>0$

$y\left(x\right)=2x+e^{-2x},x<0$

$y\left(0\right)=1,$

$y\left(1\right)=2+e^{-2},$

Exact solution: $y\left(x\right)=2x+e^{-2x}$.

Source: Cahlon and Nachman [9].

6. Discussion

The notations used in

Table 1. The comparison results for Problem 1 at $h=0.05$ and $h=0.025$.

	h
x	0.05			0.025
	2PBM5	2PBFM5	Sakai	2PBM5	2PBFM5	Sakai
0.2	1.0483$\times {10}^{-2}$	1.5847$\times {10}^{-1}$	2.1000$\times {10}^{-1}$	6.5235$\times {10}^{-4}$	6.0636$\times {10}^{-3}$	5.5000$\times {10}^{-2}$
0.4	7.6654$\times {10}^{-3}$	1.2818$\times {10}^{-1}$	1.0000$\times {10}^{-1}$	5.5278$\times {10}^{-4}$	4.8570$\times {10}^{-3}$	2.7000$\times {10}^{-2}$
0.6	6.4868$\times {10}^{-3}$	8.6763$\times {10}^{-2}$	2.0000$\times {10}^{-2}$	3.0293$\times {10}^{-4}$	3.2878$\times {10}^{-3}$	5.0000$\times {10}^{-3}$
0.8	4.3674$\times {10}^{-3}$	4.2963$\times {10}^{-2}$	2.9000$\times {10}^{-1}$	1.3336$\times {10}^{-4}$	1.6688$\times {10}^{-3}$	8.2000$\times {10}^{-2}$
1.0	2.8422$\times {10}^{-14}$	1.6875$\times {10}^{-14}$	ND	5.6843$\times {10}^{-14}$	8.8818$\times {10}^{-16}$	ND
MAXE	1.3352$\times {10}^{-2}$	1.6458$\times {10}^{-1}$	2.9000$\times {10}^{-1}$	7.2927$\times {10}^{-4}$	6.5116$\times {10}^{-3}$	8.2000$\times {10}^{-2}$
AVE	6.3798$\times {10}^{-3}$	9.1490$\times {10}^{-2}$	ND	3.6590$\times {10}^{-4}$	3.6719$\times {10}^{-3}$	ND
ITN	1	1	ND	1	1	ND
FCN	61	61	ND	101	121	ND
TS	12	12	ND	22	22	ND
Time(s)	0.000667	0.001333	ND	0.002000	0.003000	ND

,

Table 2. The comparison results for Problem 2 at $h=0.1$.

x	$\mathit{h}=0.1$
	2PBM5	2PBFM5	Bica
0.0	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
0.1	3.1147$\times {10}^{-5}$	1.93249$\times {10}^{-4}$	5.1074$\times {10}^{-5}$
0.2	3.4395$\times {10}^{-5}$	1.99280$\times {10}^{-4}$	9.5068$\times {10}^{-5}$
0.3	3.1499$\times {10}^{-5}$	1.77508$\times {10}^{-4}$	1.1406$\times {10}^{-4}$
0.4	2.6027$\times {10}^{-5}$	1.52214$\times {10}^{-4}$	1.2836$\times {10}^{-4}$
0.5	2.0952$\times {10}^{-5}$	1.17264$\times {10}^{-4}$	1.2397$\times {10}^{-4}$
0.6	1.6227$\times {10}^{-5}$	9.46777$\times {10}^{-5}$	1.1645$\times {10}^{-4}$
0.7	1.1808$\times {10}^{-5}$	6.52498$\times {10}^{-5}$	9.4795$\times {10}^{-5}$
0.8	7.6728$\times {10}^{-6}$	4.46848$\times {10}^{-5}$	7.0953$\times {10}^{-5}$
0.9	3.7449$\times {10}^{-6}$	1.93071$\times {10}^{-5}$	3.6337$\times {10}^{-5}$
1.0	1.6653$\times {10}^{-16}$	1.66533$\times {10}^{-16}$	0.0000$\times {10}^0$
MAXE	3.4395$\times {10}^{-5}$	1.9928$\times {10}^{-4}$	1.2836$\times {10}^{-4}$
AVE	1.8347$\times {10}^{-5}$	1.0634$\times {10}^{-4}$	ND
ITN	1	1	ND
FCN	29	31	ND
TS	7	7	ND
Time(s)	0.001000	0.001333	ND

,

Table 3. The comparison results for Problem 2 at $h=0.01$.

x	$\mathit{h}=0.01$
	2PBM5	2PBFM5	Bica
0.0	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
0.1	2.1026$\times {10}^{-8}$	1.6530$\times {10}^{-7}$	5.9062$\times {10}^{-7}$
0.2	1.7875$\times {10}^{-8}$	1.4016$\times {10}^{-7}$	9.4957$\times {10}^{-7}$
0.3	1.5004$\times {10}^{-8}$	1.1736$\times {10}^{-7}$	1.1637$\times {10}^{-6}$
0.4	1.2375$\times {10}^{-8}$	9.6557$\times {10}^{-8}$	1.2538$\times {10}^{-6}$
0.5	9.9553$\times {10}^{-9}$	7.7489$\times {10}^{-8}$	1.2374$\times {10}^{-6}$
0.6	7.7129$\times {10}^{-9}$	5.9896$\times {10}^{-8}$	1.1294$\times {10}^{-6}$
0.7	5.6205$\times {10}^{-9}$	4.3550$\times {10}^{-8}$	9.4221$\times {10}^{-7}$
0.8	3.6527$\times {10}^{-9}$	2.8243$\times {10}^{-8}$	6.8618$\times {10}^{-7}$
0.9	1.7863$\times {10}^{-9}$	1.3784$\times {10}^{-8}$	3.6981$\times {10}^{-7}$
1.0	4.9960$\times {10}^{-16}$	3.3307$\times {10}^{-16}$	0.0000$\times {10}^0$
MAXE	2.3421$\times {10}^{-8}$	1.8164$\times {10}^{-7}$	1.2538$\times {10}^{-6}$
AVE	1.0567$\times {10}^{-8}$	7.9674$\times {10}^{-8}$	ND
ITN	1	1	ND
FCN	209	301	ND
TS	52	52	ND
Time(s)	0.001000	0.001667	ND

,

Table 4. The comparison results for Problem 2 at $h=0.001$.

x	$\mathit{h}=0.001$
	2PBM5	2PBFM5	Bica
0.0	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
0.1	1.8201$\times {10}^{-11}$	1.5402$\times {10}^{-10}$	5.6198$\times {10}^{-9}$
0.2	1.5474$\times {10}^{-11}$	1.3057$\times {10}^{-10}$	9.4550$\times {10}^{-9}$
0.3	1.2988$\times {10}^{-11}$	1.0930$\times {10}^{-10}$	1.1765$\times {10}^{-8}$
0.4	1.0713$\times {10}^{-11}$	8.9915$\times {10}^{-11}$	1.2770$\times {10}^{-8}$
0.5	8.6178$\times {10}^{-12}$	7.2145$\times {10}^{-11}$	1.2656$\times {10}^{-8}$
0.6	6.6767$\times {10}^{-12}$	5.5756$\times {10}^{-11}$	1.1583$\times {10}^{-8}$
0.7	4.8652$\times {10}^{-12}$	4.0533$\times {10}^{-11}$	9.6804$\times {10}^{-9}$
0.8	3.1616$\times {10}^{-12}$	2.6282$\times {10}^{-11}$	7.0588$\times {10}^{-9}$
0.9	1.5459$\times {10}^{-12}$	1.2825$\times {10}^{-11}$	3.8078$\times {10}^{-9}$
1.0	4.9960$\times {10}^{-16}$	9.4369$\times {10}^{-16}$	0.0000$\times {10}^0$
MAXE	2.1113$\times {10}^{-11}$	1.7887$\times {10}^{-10}$	1.2770$\times {10}^{-8}$
AVE	9.2602$\times {10}^{-12}$	7.4977$\times {10}^{-11}$	ND
ITN	1	1	ND
FCN	2009	3001	ND
TS	502	502	ND
Time(s)	0.001667	0.006667	ND

,

Table 5. The comparison results for Problem 3 at $h=\frac{\pi }{40}$.

x	$\mathit{h}=\frac{\mathit{\pi }}{40}$
	2PBM5	2PBFM5	Bica
0.0	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
0.08	5.7369$\times {10}^{-6}$	1.5019$\times {10}^{-4}$	2.6320$\times {10}^{-5}$
0.16	5.4345$\times {10}^{-6}$	1.3806$\times {10}^{-4}$	4.3154$\times {10}^{-5}$
0.24	5.0067$\times {10}^{-6}$	1.2518$\times {10}^{-4}$	5.9168$\times {10}^{-5}$
0.31	4.4717$\times {10}^{-6}$	1.2281$\times {10}^{-4}$	6.5731$\times {10}^{-5}$
0.39	3.8703$\times {10}^{-6}$	9.3073$\times {10}^{-5}$	7.0788$\times {10}^{-5}$
0.47	3.2090$\times {10}^{-6}$	8.7944$\times {10}^{-5}$	6.6550$\times {10}^{-5}$
0.55	2.4873$\times {10}^{-6}$	5.4714$\times {10}^{-5}$	6.0163$\times {10}^{-5}$
0.63	1.7093$\times {10}^{-6}$	4.6875$\times {10}^{-5}$	4.4748$\times {10}^{-5}$
0.71	8.7848$\times {10}^{-7}$	1.0448$\times {10}^{-5}$	2.6698$\times {10}^{-5}$
0.79	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
MAXE	5.7369$\times {10}^{-6}$	1.5019$\times {10}^{-4}$	7.0788$\times {10}^{-5}$
AVE	3.2804$\times {10}^{-6}$	8.2931$\times {10}^{-5}$	ND
ITN	1	1	ND
FCN	35	31	ND
TS	7	7	ND
Time(s)	0.001000	0.001000	ND

,

Table 6. The comparison results for Problem 3 at $h=\frac{\pi }{400}$.

x	$\mathit{h}=\frac{\mathit{\pi }}{400}$
	2PBM5	2PBFM5	Bica
0.0	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
0.08	5.6732$\times {10}^{-9}$	1.6360$\times {10}^{-7}$	2.4115$\times {10}^{-7}$
0.16	5.2948$\times {10}^{-9}$	1.5256$\times {10}^{-7}$	4.3261$\times {10}^{-7}$
0.24	4.8466$\times {10}^{-9}$	1.3953$\times {10}^{-7}$	5.7246$\times {10}^{-7}$
0.31	4.3312$\times {10}^{-9}$	1.2460$\times {10}^{-7}$	6.5906$\times {10}^{-7}$
0.39	3.7513$\times {10}^{-9}$	1.0784$\times {10}^{-7}$	6.9097$\times {10}^{-7}$
0.47	3.1100$\times {10}^{-9}$	8.9341$\times {10}^{-8}$	6.6707$\times {10}^{-7}$
0.55	2.4105$\times {10}^{-9}$	6.9204$\times {10}^{-8}$	5.8648$\times {10}^{-7}$
0.63	1.6565$\times {10}^{-9}$	4.7529$\times {10}^{-8}$	4.4860$\times {10}^{-7}$
0.71	8.5168$\times {10}^{-10}$	2.4422$\times {10}^{-8}$	2.5312$\times {10}^{-7}$
0.79	4.4409$\times {10}^{-16}$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
MAXE	5.9510$\times {10}^{-9}$	1.6926$\times {10}^{-7}$	6.9097$\times {10}^{-7}$
AVE	3.4666$\times {10}^{-9}$	9.3534$\times {10}^{-8}$	ND
ITN	1	1	ND
FCN	215	301	ND
TS	52	52	ND
Time(s)	0.000667	0.001333	ND

,

Table 7. The comparison results for Problem 3 at $h=\frac{\pi }{4000}$.

x	$\mathit{h}=\frac{\mathit{\pi }}{4000}$
	2PBM5	2PBFM5	Bica
0.0	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
0.08	5.6732$\times {10}^{-12}$	1.6387$\times {10}^{-10}$	2.4100$\times {10}^{-9}$
0.16	5.2944$\times {10}^{-12}$	1.5281$\times {10}^{-10}$	4.3300$\times {10}^{-9}$
0.24	4.8457$\times {10}^{-12}$	1.3976$\times {10}^{-10}$	5.7200$\times {10}^{-9}$
0.31	4.3299$\times {10}^{-12}$	1.2480$\times {10}^{-10}$	6.5900$\times {10}^{-9}$
0.39	3.7503$\times {10}^{-12}$	1.0802$\times {10}^{-10}$	6.9100$\times {10}^{-9}$
0.47	3.1084$\times {10}^{-12}$	8.9490$\times {10}^{-11}$	6.6700$\times {10}^{-9}$
0.55	2.4079$\times {10}^{-12}$	6.9319$\times {10}^{-11}$	5.8600$\times {10}^{-9}$
0.63	1.6536$\times {10}^{-12}$	4.7608$\times {10}^{-11}$	4.4900$\times {10}^{-9}$
0.71	8.4865$\times {10}^{-13}$	2.4465$\times {10}^{-11}$	2.5300$\times {10}^{-9}$
0.79	3.5527$\times {10}^{-15}$	2.2204$\times {10}^{-15}$	0.0000$\times {10}^0$
MAXE	5.9779$\times {10}^{-12}$	1.7258$\times {10}^{-10}$	6.9100$\times {10}^{-9}$
AVE	3.4921$\times {10}^{-12}$	9.4578$\times {10}^{-11}$	ND
ITN	1	1	ND
FCN	2015	3001	ND
TS	502	502	ND
Time(s)	0.001333	0.007000	ND

and

Table 8. The comparison results for Problem 4 at $h=0.1$ and $h=0.01$.

	h
x	0.1			0.01
	2PBM5	2PBFM5	CN	2PBM5	2PBFM5	CN
0.0	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$	0.0000$\times {10}^0$
0.1	4.6295$\times {10}^{-4}$	4.6940$\times {10}^{-4}$	1.1845$\times {10}^{-3}$	6.3839$\times {10}^{-7}$	1.0880$\times {10}^{-6}$	1.4000$\times {10}^{-6}$
0.2	6.1248$\times {10}^{-4}$	6.2422$\times {10}^{-4}$	1.7120$\times {10}^{-4}$	1.0470$\times {10}^{-6}$	1.8648$\times {10}^{-6}$	2.4000$\times {10}^{-6}$
0.3	5.4674$\times {10}^{-4}$	5.6281$\times {10}^{-4}$	1.1280$\times {10}^{-4}$	1.3831$\times {10}^{-6}$	2.5046$\times {10}^{-6}$	2.9000$\times {10}^{-6}$
0.4	3.7993$\times {10}^{-4}$	4.0332$\times {10}^{-4}$	1.9790$\times {10}^{-4}$	1.6654$\times {10}^{-6}$	3.0329$\times {10}^{-6}$	3.0000$\times {10}^{-6}$
0.5	2.4926$\times {10}^{-4}$	2.6942$\times {10}^{-4}$	3.9880$\times {10}^{-4}$	1.9083$\times {10}^{-6}$	3.4671$\times {10}^{-6}$	2.8000$\times {10}^{-6}$
0.6	1.5465$\times {10}^{-4}$	1.8355$\times {10}^{-4}$	1.4130$\times {10}^{-4}$	2.1226$\times {10}^{-6}$	3.8484$\times {10}^{-6}$	2.5000$\times {10}^{-6}$
0.7	8.6837$\times {10}^{-5}$	1.0423$\times {10}^{-4}$	9.0900$\times {10}^{-5}$	2.3170$\times {10}^{-6}$	4.1361$\times {10}^{-6}$	2.0000$\times {10}^{-6}$
0.8	4.0424$\times {10}^{-5}$	6.3768$\times {10}^{-5}$	7.4100$\times {10}^{-5}$	2.4985$\times {10}^{-6}$	4.4321$\times {10}^{-6}$	1.4000$\times {10}^{-6}$
0.9	1.2431$\times {10}^{-5}$	1.4445$\times {10}^{-5}$	2.1750$\times {10}^{-4}$	2.6722$\times {10}^{-6}$	4.6167$\times {10}^{-6}$	7.0000$\times {10}^{-7}$
1.0	3.2005$\times {10}^{-7}$	2.6724$\times {10}^{-6}$	0.0000$\times {10}^0$	2.8421$\times {10}^{-6}$	4.8552$\times {10}^{-6}$	0.0000$\times {10}^0$
MAXE	6.1248$\times {10}^{-4}$	6.2422$\times {10}^{-4}$	1.1845$\times {10}^{-3}$	2.8421$\times {10}^{-6}$	4.8552$\times {10}^{-6}$	3.0000$\times {10}^{-6}$
AVE	2.5460$\times {10}^{-4}$	2.6979$\times {10}^{-4}$	ND	1.7901$\times {10}^{-6}$	3.1765$\times {10}^{-6}$	ND
ITN	2	2	ND	2	2	ND
FCN	23	31	ND	203	301	ND
TS	7	7	ND	52	52	ND
Time(s)	0.000333	0.001000	ND	0.001333	0.001667	ND

can be found in “Abbreviations” at the end of this paper.

Table 1, Table 2, Table 3, Table 4, Table 5, Table 6 and Table 7 demonstrate that 2PBM5 gives more accurate results compared to 2PBFM5, Sakai and Bica for both constant and pantograph delay type problems. We can observe that the absolute errors for 2PBM5 at each point x is smaller than 2PBFM5 [14], Sakai [5], and Bica [10] as the step size, h decreases. On the other hand, Table 8 shows that 2PBM5 is in a good agreement with previous methods, 2PBFM5 and CN, in solving time-dependent delay type problem, where 2PBM5 slightly improved in accuracy when $h=0.1$ and $h=0.01$ for maximum absolute error (MAXE).

The total function calls in Table 1, Table 2, Table 3, Table 4, Table 5, Table 6, Table 7 and Table 8 for 2PBM5 are smaller than 2PBFM5 in most of the tested problems. This is due to the 2PBM5 being in a diagonally implicit approach where the function calls needed for every iteration are slightly smaller than the fully implicit approach in 2PBFM5. We also observed that the execution time for 2PBM5 is faster than 2PBFM5 in Table 1, Table 2, Table 3, Table 4, Table 5, Table 6, Table 7 and Table 8.

The total steps and number of guesses for both multistep methods, 2PBM5 and 2PBFM5, are the same since both are block methods that calculate two solutions simultaneously at every iteration and solved using the same strategy of the Newton-like method in shooting technique.

Figure 3 displays the performance of the 2PBM5 method compared to existing methods subject to accuracy as the h decreased.

7. Conclusions

In this paper, we have shown that the proposed method 2PBM5 can produce better accuracy compared to 2PBFM5 [14], Sakai [5], Bica [10], and CN [9] for solving constant and non-constant delay types. The 2PBM5 is also faster in computational time compared to 2PBFM5. The advantages of the diagonal form for the formula in 2PBM5 reflect the total function calls where it is lesser than 2PBFM5. However, the limitation of the proposed method in this research is only for solving the DDEs of non-discontinuity in solution cases. In conclusion, the numerical simulation has shown that 2PBM5 is computationally efficient in solving the second-order DDEs of the constant, pantograph, and time-dependent delay type with boundary conditions.