Numerical Solution for Singular Boundary Value Problems Using a Pair of Hybrid Nyström Techniques

Abstract

This manuscript presents an efficient pair of hybrid Nyström techniques to solve second-order Lane–Emden singular boundary value problems directly. One of the proposed strategies uses three off-step points. The obtained formulas are paired with an appropriate set of formulas implemented for the first step to avoid singularity at the left end of the integration interval. The fundamental properties of the proposed scheme are analyzed. Some test problems, including chemical kinetics and physical model problems, are solved numerically to determine the efficiency and validity of the proposed approach.

1. Introduction

In this paper, we consider the following two-point singular boundary value problem (SBVP):

$$y^{″}\left(x\right)+\frac{\lambda }{x}{y}^{\prime }\left(x\right)=k(x,y),0<x\le x_N=1.$$

(1)

We consider this together with any of the two-point boundary conditions:

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

(2)

or

$$y\left(0\right)=y_a,y^{\prime }\left(1\right)=y_b^{\prime },$$

(3)

or

$$y^{\prime }\left(0\right)=y_a^{\prime },y\left(1\right)=y_b,$$

(4)

where $\lambda ,y_a,y_b,y_a^{\prime },y_b^{\prime },$ are known real values and $k(x,y)$ denotes a continuous real function, where we assume the necessary conditions to guarantee the existence of a unique solution to the problem. The existence and uniqueness of the solution to the problem (1) subjected to any of the boundary conditions above have been established by Pandey [1] and Zou [2].

According to Thula and Roul [3], the mathematical expression of numerous problems arising in chemical kinetics, astrophysics, catalytic diffusion reactions, celestial mechanics, engineering, and various physical models gives rise to second-order singular boundary value ordinary differential equations (SSBODEs) of the type given in (1).

The problem of reactant concentration in a chemical reactor, reaction–diffusion processes inside a porous catalyst, the conduction of heat in the human head, the distribution of oxygen in a spherical shell, and many others, can be modelled by the system (1).

Significantly, much work has been carried out to obtain numerical solutions for the above singular problems. Different strategies have been reported for solving (1), where the fundamental difficulty arises due to the singularity at $x=0$. Notable scholars in the field of numerical analysis have proposed numerical techniques for solving the problem (1). Examples of such techniques include the finite difference methods (FDM) proposed in [4,5], the spline methods (SM) proposed in [6,7], the Padè approximation method (PAM) introduced in [8,9], the pseudospectral method (PM) proposed in [10], or the Jacobi–Gauss collocation method (JCM) reported in [11]. Other manuscripts on recently developed numerical or analytical techniques for solving (1)–(4) inlcude those of [12,13,14,15,16,17,18].

We propose a pair of hybrid Nyström techniques (PHNT) for solving the SBVPs given in (1) numerically. The main formulas provide a method of order six, which cannot solve the problem on its own due to the singularity at $x=0$. We have designed a second ad hoc method that applies only to the first subinterval and is unaffected by singularity. In this way, we obtain a scheme capable of solving the problem posed effectively. Comparisons show that the proposed method is more advantageous than other existing methods. A problem of particular interest is how to select the optimal value of h. We have not addressed this issue here, but the CESTAC method and the CADNA library could be helpful for this task [19,20]. This technique is based on the use of stochastic arithmetic in place of floating-point arithmetic to validate the results and find an optimal solution.

The present work is outlined as follows. In Section 2, we present the PHNT method for solving SBVPs. The characteristics of the developed formulas are analyzed in Section 3. Some issues with the implementation of the PHNT are considered in Section 4. In Section 5, we present the numerical results of some physical models and catalytic diffusion–reaction problems to show the efficiency and reliability of the proposed technique. Some conclusions are outlined in Section 6.

2. Development of the PHNT Method

To obtain the PHNT method, we firstly reformulate the equation in (1) as $y^{″}\left(x\right)=f(x,y\left(x\right),y^{\prime }\left(x\right))$, where $f(x,y\left(x\right),y^{\prime }\left(x\right))=k(x,y\left(x\right))-\frac{\lambda }{x}{y}^{\prime }\left(x\right)$. Thus, the singularity is transferred to the function f.

2.1. Main Formulas

We consider that the exact solution $y\left(x\right)$ of the SBVP on an interval $[x_n,x_{n+1}]$, $x_n>0$ with step size $h=x_{n+1}-x_n$ can be approximated by a polynomial $p\left(x\right)$ in the form:

$$y\left(x\right)\simeq p\left(x\right)=\sum _{j=0}^6{a}_j{x}^j,$$

(5)

From this, it readily follows that:

$$y^{\prime }\left(x\right)\simeq p^{\prime }\left(x\right)=\sum _{j=1}^6{a}_{j}j{x}^{j-1},$$

(6)

$$y^{″}\left(x\right)\simeq p^{″}\left(x\right)=\sum _{j=2}^6{a}_{j}j(j-1)x^{j-2},$$

(7)

where $a_j\in \mathbb{R}$ are unknown coefficients that must be specified using collocation conditions at some chosen nodes.

We take the following intermediate nodes on $[x_n,x_{n+1}]$: $x_{n+u}=x_n+uh,x_{n+v}=x_n+vh$ and $x_{n+w}=x_n+wh$ with $0<u<v<w<1$. Consider the approximations in (5) and (6) evaluated at $x_n$, as well as the one in (7) evaluated at $x_n,x_{n+u},x_{n+v},x_{n+w},x_{n+1}$. Doing this, we obtain a system of equations with seven unknowns $a_n,n=0\left(1\right)6$, given by:

$$\begin{array}{ccc}& & p\left(x_n\right)=y_n,p^{\prime }\left(x_n\right)=y_n^{\prime },p^{″}\left(x_n\right)=f_n,\hfill \\ & & p^{″}\left(x_{n+u}\right)=f_{n+u},p^{″}\left(x_{n+v}\right)=f_{n+v},p^{″}\left(x_{n+w}\right)=f_{n+w},p^{″}\left(x_{n+1}\right)=f_{n+1},\hfill \end{array}$$

where $y_{n+j}$ and $f_{n+j}$ denote approximations of $y\left(x_{n+j}\right)$ and $y^{″}\left(x_{n+j}\right)$, respectively. The system may be written in matrix form as:

$$\left(\begin{array}{ccccccc}1& x_n& x_n^2& x_n^3& x_n^4& x_n^5& x_n^6\\ 0& 1& 2x_n& 3x_n^2& 4x_n^3& 5x_n^4& 6x_n^5\\ 0& 0& 2& 6x_n& 12x_n^2& 20x_n^3& 30x_n^4\\ 0& 0& 2& 6x_{n+u}& 12x_{n+u}^2& 20x_{n+u}^3& 30x_{n+u}^4\\ 0& 0& 2& 6x_{n+v}& 12x_{n+v}^2& 20x_{n+v}^3& 30x_{n+v}^4\\ 0& 0& 2& 6x_{n+w}& 12x_{n+w}^2& 20x_{n+w}^3& 30x_{n+w}^4\\ 0& 0& 2& 6x_{n+1}& 12x_{n+1}^2& 20x_{n+1}^3& 30x_{n+1}^4\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6\end{array}\right)=\left(\begin{array}{c}{y}_n\\ y_n^{\prime }\\ f_n\\ f_{n+u}\\ f_{n+v}\\ f_{n+w}\\ f_{n+1}\end{array}\right).$$

Solving this system, we obtain the values of $a_n,n=0\left(1\right)6$. Using the substitution $x=x_n+zh,$ the polynomial in (5) may be expressed as:

$$\begin{array}{ccc}\hfill p(x_n+zh)& =& {\alpha }_0\left(z\right)y_n+h{\alpha }_1\left(z\right)y_n^{\prime }\hfill \\ & & +h^2({\beta }_0\left(z\right)f_n+{\beta }_u\left(z\right)f_{n+u}+{\beta }_v\left(z\right)f_{n+v}+{\beta }_w\left(z\right)f_{n+w}+{\beta }_1\left(z\right)f_{n+1}),\hfill \end{array}$$

(8)

where $u=\frac{1}{3},v=\frac{1}{2},w=\frac{2}{3}$, and the coefficients ${\alpha }_0\left(z\right)=1$, ${\alpha }_1\left(z\right)$, and ${\left\{{\beta }_i\left(z\right)\right\}}_{i=0,u,v,w,1}$ depend on z.

Evaluating the formula in (8) and its derivative at $z=1$, we obtain approximations of $y\left(x_{n+1}\right)$ and $y^{\prime }\left(x_{n+1}\right)$, given, respectively, by:

$$\begin{array}{ccc}& & y_{n+1}=y_n+h{y}_n^{\prime }+h^2\left(\frac{9f_{n+u}}{20}-\frac{4f_{n+v}}{15}+\frac{9f_{n+w}}{40}+\frac{11f_n}{120}\right),\hfill \\ & & y_{n+1}^{\prime }=y_n^{\prime }+\frac{h}{120}\left(11f_n+81f_{n+u}+81f_{n+v}-64f_{n+w}+11f_{n+1}\right).\hfill \end{array}$$

(9)

Now, evaluating $p\left(x\right)$ and $p^{\prime }\left(x\right)$ at $x_{n+u},x_{n+v},x_{n+w}$, we obtain the following hybrid Nyström-type formulas:

$$\begin{array}{ccc}\hfill y_{n+u}& =& y_n-\frac{1}{3}h{y}_n^{\prime }+\frac{h^2\left(240f_{n+u}-224f_{n+v}+87f_{n+v}+83f_n-6f_{n+1}\right)}{3240},\hfill \\ \hfill y_{n+v}& =& y_n+\frac{h}{2}{y}_n^{\prime }+h^2\left(\frac{207f_{n+u}}{1280}-\frac{f_{n+v}}{8}+\frac{63f_{n+w}}{1280}+\frac{163f_n}{3840}-\frac{13f_{n+1}}{3840}\right),\hfill \\ \hfill y_{n+w}& =& y_n+\frac{2h}{3}{y}_n^{\prime }+\frac{2h^2}{405}\left(51f_{n+u}-32f_{n+v}+15f_{n+w}+12f_n-f_{n+1}\right).\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill y_{n+u}^{\prime }& =& y_n^{\prime }+h\left(\frac{19f_{n+u}}{40}-\frac{152f_{n+v}}{405}+\frac{17f_{n+w}}{120}+\frac{329f_n}{3240}-\frac{31f_{n+1}}{3240}\right),\hfill \\ \hfill y_{n+v}^{\prime }& =& y_n^{\prime }+\frac{h\left(1053f_{n+u}-512f_{n+v}+243f_{n+w}+193f_n-17f_{n+1}\right)}{1920},\hfill \\ \hfill y_{n+w}^{\prime }& =& y_n^{\prime }+\frac{1}{405}h\left(216f_{n+u}-64f_{n+v}+81f_{n+w}+41f_n-4f_{n+1}\right).\hfill \end{array}$$

(11)

The formulas in (9)–(11) altogether form the main block method.

2.2. Formulas to Circumvent the Singularity

This block method cannot be used directly for solving a BVP problem with the differential equation in $\left(1\right)$ because it is not possible to evaluate $f_0=f(x_0,y_0,y_0^{\prime })$, since there is a singularity at $x_0=0$.

To overcome this drawback, we have developed a set of formulas specially designed for the subinterval $[x_0,x_1]$, where the value $f_0$ is absent. These formulas are obtained similarly to before as:

$$\begin{array}{ccc}\hfill y_1& =& y_0+h{y}_0^{\prime }+h^2\left(\frac{51f_u}{40}-\frac{26f_v}{15}+\frac{21f_w}{20}-\frac{1}{120}11f_1\right),\hfill \\ \hfill y_1^{\prime }& =& y_0^{\prime }+h\left(\frac{13f_u}{9}-\frac{16f_v}{9}+\frac{10f_w}{9}-\frac{f_1}{9}\right).\hfill \end{array}$$

(12)

For the remaining formulas, we obtain:

$$\begin{array}{ccc}\hfill y_u& =& y_0+\frac{h}{3}{y}_0^{\prime }+h^2\left(\frac{329f_u}{1080}-\frac{194f_v}{405}+\frac{139f_w}{540}-\frac{89f_1}{3240}\right),\hfill \\ \hfill y_v& =& y_0+\frac{h}{2}{y}_0^{\prime }+h^2\left(\frac{87f_u}{160}-\frac{193f_v}{240}+\frac{69f_w}{160}-\frac{1}{240}11f_1\right),\hfill \\ \hfill y_w& =& y_0+\frac{2h}{3}{y}_0^{\prime }-\frac{1}{405}2h^2\left(-159f_u+224f_v-123f_w+13f_1\right).\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill y_u^{\prime }& =& y_0^{\prime }-\frac{1}{18}h\left(-25f_u+36f_v-19f_w+2f_1\right),\hfill \\ \hfill y_v^{\prime }& =& y_0^{\prime }-\frac{1}{64}h\left(-93f_u+120f_v-66f_w+7f_1\right),\hfill \\ \hfill y_w^{\prime }& =& y_0^{\prime }-\frac{1}{9}h\left(-13f_u+16f_v-10f_w+f_1\right).\hfill \end{array}$$

(14)

Taking a small step size h, and considering all the formulas in (9)–(11) for $n=1,2,\dots ,N-1$, together with the ones developed in (12)–(14) for the first step, we obtain a global method that can provide accurate approximations to complete the integration along the interval $[0,x_N]$.

3. Characteristics of the Method

The main properties of the proposed technique PHNT are studied here, where the most challenging task is to analyze the convergence of the global method.

3.1. Consistency and Order of the Formulas

The formulas in $\left(9\right)$–$\left(11\right)$ may be written as:

$$\begin{array}{c}\hfill \overline{A}{V}_n=h\overline{B}{V}_n^{\prime }+h^2\overline{D}{F}_n,\end{array}$$

(15)

where $\overline{A},\overline{B},\overline{D}$ are constant matrices containing the coefficients of the formulas (9)–(11), and:

$$\begin{array}{ccc}& & V_n={\left(y_n,y_{n+u},y_{n+v},y_{n+w},y_{n+1}\right)}^T,\hfill \\ & & V_n^{\prime }={\left(y_n^{\prime },y_{n+u}^{\prime },y_{n+v}^{\prime },y_{n+w}^{\prime },y_{n+1}^{\prime }\right)}^T,\hfill \\ & & F_n={\left(f_n,f_{n+u},f_{n+v},f_{n+w},f_{n+1}\right)}^T.\hfill \end{array}$$

Using standard strategies (see [21]), assuming that $y\left(x\right)$ has enough derivatives, we define the operator ℓ related to the formulas in $\left(9\right)$–$\left(11\right)$:

$$\begin{array}{ccc}\hfill \ell \left[y\right(x);h]& =& \sum _{j\in I}\left[{\alpha }_{j}y\left(x_n+jh\right)-h{\beta }_j{y}^{\prime }\left(x_n+jh\right)-h^2{\gamma }_j{y}^{″}\left(x_n+jh\right)\right],\hfill \end{array}$$

(16)

where ${\alpha }_j,{\beta }_j$, and ${\gamma }_j$ are, respectively, vector columns of $\overline{A}, \overline{B}$ and $\overline{D}$, while I denotes the set of indices, $I=\{0,u,v,w,1\}$. Expanding in Taylor series about $x_n$, we obtain:

$$\ell [y\left(x\right);h]=C_{0}y\left(x_n\right)+C_{1}h{y}^{\prime }\left(x_n\right)+C_2{h}^2{y}^{″}\left(x_n\right)+\dots +C_q{h}^q{y}^q\left(x_n\right)+\dots ,$$

(17)

where:

$$\begin{array}{ccc}\hfill C_q& =& \frac{1}{q!}\left[\sum _{j\in I}^k{j}^q{\alpha }_j-q\sum _{j\in I}^k{j}^{q-1}{\beta }_j-q(q-1)\sum _{j\in I}^k{j}^{q-2}{\gamma }_j\right],q=0,1,2,3,\dots .\hfill \end{array}$$

(18)

According to [16], the above operator and the associated formulas are said to be of order p if $C_0=C_1=...=C_{p+1}=0,C_{p+2}\ne 0$, with $C_{p+2}$ as the vector of local truncation errors. For the formulas in (9)–(11), we obtain $C_0=C_1=\dots =C_6=0$ and:

$$C_7={\left(\frac{1}{\mathrm{181,440}},-\frac{1}{\mathrm{1,088,640}},\frac{67}{\mathrm{44,089,920}},\frac{1}{\mathrm{362,880}},\frac{11}{\mathrm{2,755,620}},\frac{1}{\mathrm{131,220}},\frac{1}{\mathrm{138,240}},\frac{1}{\mathrm{131,220}}\right)}^T,$$

This shows that each of the above formulas is of order 5. Since the order of the formulas is greater than one, they are consistent. For the ad hoc formulas used for the first step, it is easy to see that they are also consistent.

3.2. Convergence Analysis

We start by defining convergence, then we will show that the proposed method is convergent by writing all the formulas in (9)–(14) in an appropriate matrix-vector form.

Definition 1.

Let $y\left(x\right)$ denote the exact solution of the given singular boundary value problem and let ${\left\{y_j\right\}}_{j=0}^N$ be the approximations obtained with the developed numerical strategy. The method is said to be convergent of order p if, for a sufficiently small h, there exists a constant C independent of h, such that:

$$\underset{0\le j\le N}{max}|y\left(x_j\right)-y_j|\le C{h}^p.$$

Note that in this situation, we obtain ${\displaystyle \underset{0\le j\le N}{max}|y\left(x_j\right)-y_j|\to 0}$ as $h\to 0$.

Theorem 1

(Convergence theorem). Let $y\left(x\right)$ be the true solution of the SBVP in $\left(1\right)$ with the boundary conditions in $\left(2\right)$, and ${\left\{y_j\right\}}_{j=0}^N$ the discrete solution provided by the proposed global method. Then, the proposed method is convergent to order six.

Proof. 

Following [22], we define the matrix D of dimension $8N\times 8N$ given by:

$$D=\left[\begin{array}{cccc}{D}_{1,1}& D_{1,2}& \dots & D_{1,2N}\\ ⋮& ⋮& & ⋮\\ D_{2N,1}& D_{2N,2}& \dots & D_{2N,2N}\end{array}\right],$$

where the elements $D_{i,j}$ are $4\times 4$ sub-matrices, except the $D_{i,N+1},i=1,\dots ,2N$, which have a size of $4\times 3$, and $D_{i,2N},i=1,\dots ,2N$, which have a size of $4\times 5$. Those sub-matrices are

$$D_{i,i}=I,i=N+2,\dots ,2N,\mathrm{being}I\mathrm{the}\mathrm{identity}\mathrm{matrix}\mathrm{of}\mathrm{order}\mathrm{four},$$

$$D_{N,N}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right];D_{i,i-1}=\left[\begin{array}{cccc}0& 0& 0& -1\\ 0& 0& 0& -1\\ 0& 0& 0& -1\\ 0& 0& 0& -1\end{array}\right],i=N+2,\dots ,2N;$$

$$D_{N+1,N+1}=\left[\begin{array}{ccccc}-1& 1& 0& 0& 0\\ -1& 0& 1& 0& 0\\ -1& 0& 0& 1& 0\\ -1& 0& 0& 0& 1\end{array}\right];D_{1,N+1}=h\left[\begin{array}{ccccc}-{\displaystyle \frac{1}{3}}& 0& 0& 0& 0\\ -{\displaystyle \frac{1}{2}}& 0& 0& 0& 0\\ -{\displaystyle \frac{2}{3}}& 0& 0& 0& 0\\ -1& 0& 0& 0& 0\end{array}\right];$$

$$D_{i,N+i}=h\left[\begin{array}{cccc}0& 0& 0& -{\displaystyle \frac{1}{3}}\\ 0& 0& 0& -{\displaystyle \frac{1}{2}}\\ 0& 0& 0& -{\displaystyle \frac{2}{3}}\\ 0& 0& 0& -1\end{array}\right],i=1\dots ,N-1;D_{N,2N}=h\left[\begin{array}{ccccc}0& 0& 0& 0& -{\displaystyle \frac{1}{3}}\\ 0& 0& 0& 0& -{\displaystyle \frac{1}{2}}\\ 0& 0& 0& 0& -{\displaystyle \frac{2}{3}}\\ 0& 0& 0& 0& -1\end{array}\right].$$

For the rest of the submatrices not included above, it is $D_{i,j}=\mathbb{O}$—that is, they are null matrices.

We also define the matrix U of dimension $8N\times (4N+1)$:

$$U=\left[\begin{array}{cccc}{U}_{1,1}& U_{1,2}& \dots & U_{1,N}\\ ⋮& ⋮& & ⋮\\ U_{2N,1}& U_{2N,2}& \dots & U_{2N,N}\end{array}\right],$$

where the elements $U_{i,j}$ are $4\times 4$ submatrices except the $U_{i,1},i=1,\dots ,2N,$ which have a size of $4\times 5$. These submatrices are given as follows:

$$U_{1,1}=h\left[\begin{array}{ccccc}0& -\frac{329}{1080}& \frac{194}{405}& -\frac{139}{540}& \frac{89}{3240}\\ 0& -\frac{87}{160}& \frac{193}{240}& -\frac{69}{160}& \frac{11}{240}\\ 0& -\frac{106}{135}& \frac{448}{405}& -\frac{82}{135}& \frac{26}{405}\\ 0& -\frac{51}{40}& \frac{26}{15}& -\frac{21}{20}& \frac{11}{120}\end{array}\right];$$

$$U_{i,i}=h\left[\begin{array}{cccc}-\frac{2}{27}& \frac{28}{405}& -\frac{29}{1080}& \frac{1}{540}\\ -\frac{207}{1280}& \frac{1}{8}& -\frac{63}{1280}& \frac{13}{3840}\\ -\frac{34}{135}& \frac{64}{405}& -\frac{2}{27}& \frac{2}{405}\\ -\frac{9}{20}& \frac{4}{15}& -\frac{9}{40}& 0\end{array}\right],i=2\dots ,N;$$

$$U_{i,i-1}=h\left[\begin{array}{cccc}0& 0& 0& -\frac{83}{3240}\\ 0& 0& 0& -\frac{163}{3840}\\ 0& 0& 0& -\frac{8}{135}\\ 0& 0& 0& -\frac{11}{120}\end{array}\right],i=3,\dots ,N;U_{2,1}=h\left[\begin{array}{ccccc}0& 0& 0& 0& -\frac{83}{3240}\\ 0& 0& 0& 0& -\frac{163}{3840}\\ 0& 0& 0& 0& -\frac{8}{135}\\ 0& 0& 0& 0& -\frac{11}{120}\end{array}\right];$$

$$U_{N+1,1}=\left[\begin{array}{ccccc}0& -\frac{25}{18}& 2& -\frac{19}{18}& \frac{1}{9}\\ 0& -\frac{93}{64}& \frac{15}{8}& -\frac{33}{32}& \frac{7}{64}\\ 0& -\frac{13}{9}& \frac{16}{9}& -\frac{10}{9}& \frac{1}{9}\\ 0& -\frac{3}{2}& 2& -\frac{3}{2}& 0\end{array}\right];$$

$$U_{N+j,j}=\left[\begin{array}{cccc}-\frac{19}{40}& \frac{152}{405}& -\frac{17}{120}& \frac{31}{3240}\\ -\frac{351}{640}& \frac{4}{15}& -\frac{81}{640}& \frac{17}{1920}\\ -\frac{8}{15}& \frac{64}{405}& -\frac{1}{5}& \frac{4}{405}\\ -\frac{27}{40}& \frac{8}{15}& -\frac{27}{40}& -\frac{11}{120}\end{array}\right],j=2,\dots ,N;$$

$$U_{N+j,j-1}=\left[\begin{array}{cccc}0& 0& 0& -\frac{329}{3240}\\ 0& 0& 0& -\frac{193}{1920}\\ 0& 0& 0& -\frac{41}{405}\\ 0& 0& 0& -\frac{11}{120}\end{array}\right],j=3,\dots ,N;U_{N+2,1}=\left[\begin{array}{ccccc}0& 0& 0& 0& -\frac{329}{3240}\\ 0& 0& 0& 0& -\frac{193}{1920}\\ 0& 0& 0& 0& -\frac{41}{405}\\ 0& 0& 0& 0& -\frac{11}{120}\end{array}\right].$$

For the rest of submatrices $U_{i,j}$ not included above, it is $U_{i,j}=\mathbb{O}$—that is, they are null matrices.

We note that the submatrices $D_{i,j}$ and $U_{i,j}$ contain the coefficients of the formulas in (12)–(14) and those of the formulas in (9)–(11), for $n=1,2,\dots ,N-1$.

Let us denote the vectors of exact values as:

$$\begin{array}{ccc}\hfill Y& =& {\left(y\left(x_u\right),y\left(x_v\right),y\left(x_w\right),y\left(x_1\right),\dots ,y\left(x_{N-1+w}\right),y^{\prime }\left(x_0\right),y^{\prime }\left(x_u\right),\dots ,y^{\prime }\left(x_N\right)\right)}^T,\hfill \\ \hfill F& =& \left(f(x_0,y\left(x_0\right),y^{\prime }\left(x_0\right)),f(x_u,y\left(x_u\right),y^{\prime }\left(x_u\right)),\dots ,f(x_N,y\left(x_N\right),y^{\prime }\left(x_N\right)))\right..\hfill \end{array}$$

Note that Y has $(4N-1)+(4N+1)=8N$ components while F has $(4N+1)$ components, because, due to the boundary conditions in (2), $y\left(x_0\right)$ and $y\left(x_N\right)$ are known values, $y\left(x_0\right)=y_a,y\left(x_N\right)=y_b$.

The exact form of the discretized formulas to approximate the boundary value problem can be written as:

$$D_{8N\times 8N}{Y}_{8N}+h{U}_{8N\times (4N+1)}{F}_{4N+1}+C_{8N}=L{\left(h\right)}_{8N},$$

(19)

where we have included the dimensions for clarity. Here, $C_{8N}$ is a vector that contains the known values—that is:

$$C_{8N}={\left(-y_a,-y_a,-y_a,-y_a,0,\dots ,0,y_b,0,\dots ,0\right)}^T,$$

and $L{\left(h\right)}_{8N}$ is another vector containing the LTEs of the formulas, given by:

$$L{\left(h\right)}_{8N}\simeq \left[\begin{array}{c}\frac{83h^6{y}^{\left(6\right)}\left(x_0\right)}{\mathrm{699,840}}+\mathcal{O}\left(h^7\right)\\ \frac{163h^6{y}^{\left(6\right)}\left(x_0\right)}{\mathrm{829,440}}+\mathcal{O}\left(h^7\right)\\ \frac{h^6{y}^{\left(6\right)}\left(x_0\right)}{3645}+\mathcal{O}\left(h^7\right)\\ \frac{11h^6{y}^{\left(6\right)}\left(x_0\right)}{\mathrm{25,920}}+\mathcal{O}\left(h^9\right)\\ \frac{67h^7{y}^{\left(7\right)}\left(x_1\right)}{\mathrm{44,089,920}}+\mathcal{O}\left(h^8\right)\\ \frac{h^7{y}^{\left(7\right)}\left(x_1\right)}{\mathrm{362,880}}+\mathcal{O}\left(h^9\right)\\ \frac{11h^7{y}^{\left(7\right)}\left(x_1\right)}{\mathrm{2,755,620}}+\mathcal{O}\left(h^8\right)\\ \frac{h^7{y}^{\left(7\right)}\left(x_1\right)}{\mathrm{181,440}}+\mathcal{O}\left(h^{10}\right)\\ ⋮\\ \frac{h^7{y}^{\left(7\right)}\left(x_{N-1}\right)}{\mathrm{181,440}}+\mathcal{O}\left(h^8\right)\\ \frac{329h^5{y}^{\left(6\right)}\left(x_0\right)}{\mathrm{699,840}}+\mathcal{O}\left(h^7\right)\\ \frac{193h^5{y}^{\left(6\right)}\left(x_0\right)}{\mathrm{414,720}}+\mathcal{O}\left(h^7\right)\\ \frac{41h^5{y}^{\left(6\right)}\left(x_0\right)}{\mathrm{87,480}}+\mathcal{O}\left(h^7\right)\\ \frac{11h^5{y}^{\left(6\right)}\left(x_0\right)}{\mathrm{25,920}}+\mathcal{O}\left(h^7\right)\\ \frac{h^6{y}^{\left(7\right)}\left(x_1\right)}{\mathrm{131,220}}+\mathcal{O}\left(h^7\right)\\ \frac{h^6{y}^{\left(7\right)}\left(x_1\right)}{\mathrm{138,240}}+\mathcal{O}\left(h^7\right)\\ \frac{h^6{y}^{\left(7\right)}\left(x_1\right)}{\mathrm{131,220}}+\mathcal{O}\left(h^7\right)\\ -\frac{h^7{y}^{\left(8\right)}\left(x_1\right)}{\mathrm{1,088,640}}+\mathcal{O}\left(h^8\right)\\ ⋮\\ -\frac{h^{10}{y}^{\left(8\right)}\left(x_{N-1}\right)}{\mathrm{1,088,640}}+\mathcal{O}\left(h^{11}\right)\end{array}\right].$$

Concerning the approximate values, they are provided by the system:

$$D_{8N\times 8N}{\overline{Y}}_{8N}+h{U}_{8N\times (4N+1)}{\overline{F}}_{4N+1}+C_{8N}=0,$$

(20)

where ${\overline{Y}}_{8N}$ approximates the vector $Y_{8N}$—that is:

$${\overline{Y}}_{8N}={\left(y_u,y_v,y_w,y_1,\dots ,y_{N-1+w},y_0^{\prime },y_u^{\prime },\dots ,y_N^{\prime }\right)}^T,$$

In addition:

$${\overline{F}}_{4N+1}={\left(f_0,f_u,f_v,f_w,f_1,\dots ,f_N\right)}^T.$$

We subtract (20) from (19) to obtain:

$$D_{8N\times 8N}{E}_{8N}+h{U}_{8N\times (4N+1)}{\left(F-\overline{F}\right)}_{4N+1}=L{\left(h\right)}_{8N},$$

(21)

where $E_{8N}=Y_{8N}-{\overline{Y}}_{8N}={\left(e_u,e_v,\dots ,e_{N-1+w},e_0^{\prime },e_u^{\prime },\dots ,e_N^{\prime }\right)}^T$ is a vector of errors at the discrete points.

Through the Mean-Value Theorem in [23], we can put any convenient subindex i as:

$$\begin{array}{ccc}& & f(x_i,y\left(x_i\right),y^{\prime }\left(x_i\right))-f(x_i,y_i,y_i^{\prime })=\left(y\left(x_i\right)-y_i\right)\frac{\partial f}{\partial y}\left({\xi }_i\right)+\left(y^{\prime }\left(x_i\right)-y_i^{\prime }\right)\frac{\partial f}{\partial y^{\prime }}\left({\xi }_i\right),\hfill \end{array}$$

where ${\xi }_i$ denotes intermediate points in the line between $(x_i,y\left(x_i\right),y^{\prime }\left(x_i\right))$ and $(x_i,y_i,y_i^{\prime })$. Thus, we can say that:

$$\begin{array}{ccc}\hfill {(F-\overline{F})}_{4N+1}& =& \left(\begin{array}{cccccccc}\frac{\partial f}{\partial y}\left({\xi }_0\right)& 0& \dots & 0& \frac{\partial f}{\partial y^{\prime }}\left({\xi }_0\right)& 0& \dots & 0\\ 0& \frac{\partial f}{\partial y}\left({\xi }_u\right)& \dots & 0& 0& \frac{\partial f}{\partial y^{\prime }}\left({\xi }_u\right)& \dots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \frac{\partial f}{\partial y}\left({\xi }_N\right)& 0& 0& \dots & \frac{\partial f}{\partial y^{\prime }}\left({\xi }_N\right)\end{array}\right)\left(\begin{array}{c}{e}_0\\ e_u\\ ⋮\\ e_N\\ e_0^{\prime }\\ e_u^{\prime }\\ ⋮\\ e_N^{\prime }\end{array}\right)\hfill \\ & =& \left(\begin{array}{cccccccc}0& \dots & 0& \frac{\partial f}{\partial y^{\prime }}\left({\xi }_0\right)& 0& \dots & 0& 0\\ \frac{\partial f}{\partial y}\left({\xi }_u\right)& \dots & 0& 0& \frac{\partial f}{\partial y^{\prime }}\left({\xi }_u\right)& \dots & 0& 0\\ ⋮& \ddots & 0& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& \dots & \frac{\partial f}{\partial y}\left({\xi }_{N-1+w}\right)& 0& 0& \dots & \frac{\partial f}{\partial y^{\prime }}\left({\xi }_{N-1+w}\right)& 0\\ 0& \dots & 0& 0& 0& \dots & 0& \frac{\partial f}{\partial y^{\prime }}\left({\xi }_N\right)\end{array}\right)\left(\begin{array}{c}{e}_u\\ e_v\\ ⋮\\ e_{N-1+w}\\ e_0^{\prime }\\ e_u^{\prime }\\ ⋮\\ e_N^{\prime }\end{array}\right)\hfill \\ & =& J_{(4N+1)\times 8N}{E}_{8N}.\hfill \end{array}$$

Note that in the second identity we have used, $e_0=y\left(x_0\right)-y_a=0$ and $e_N=y\left(x_N\right)-y_b=0$.

In view of the above, the equation in (21) may be arranged as:

$$\left(D_{8N\times 8N}+h{U}_{8N\times (4N+1)}{J}_{(4N+1)\times 8N}\right)E_{8N}=L{\left(h\right)}_{8N},$$

(22)

Setting $M=D+hUJ$, we simply find that:

$$M_{8N\times 8N}{E}_{8N}=L{\left(h\right)}_{8N}.$$

(23)

Following [24], we prove that, except for a few selected values of $h>0$, matrix M is invertible. If we use the abbreviated notation $D_N=D_{8N\times 8N}$, given the form of this matrix where the submatrices have many zeros, it is easy to verify that, for $N=2$, the determinant is $|D_2|=-2h$. Now, by induction, it can be proven that $|D_N|=-Nh$; thus $D_N$ is invertible as long as it is $h>0$.

Now, the matrix M may be rewritten as:

$$M=D+hUJ=(Id-B)D,$$

where $Id$ is the identity matrix of order $8N$ and $B=-hUJ{D}^{-1}$. Thus, we find that $\left|M\right|=|Id-B|\left|D\right|$.

As $|\lambda Id-B|={\prod }_{i=1}^{8N}(\lambda -{\lambda }_i)$ is the characteristic polynomial of B, in order to have $|Id-B|\ne 0$, if we take $\lambda =1$ it is sufficient to choose h, such that:

$$\left\{1/{\overline{\lambda }}_i:{\overline{\lambda }}_i \mathrm{is\; an\; eigenvalue\; of} UJ{D}^{-1}\right\}.$$

For such values of h, the equation in (23) may be rewritten as:

$$E=\left(M^{-1}\right)L\left(h\right).$$

(24)

The maximum norm in $\mathbb{R}$, ${\displaystyle \parallel E\parallel =\underset{i}{max}\left|e_i\right|}$, and the corresponding matrix-induced norm in ${\mathbb{R}}^{8N\times 8N}$ are considered. If we expand the terms of $M^{-1}$ in powers of h, it can be shown that $\parallel M^{-1}\parallel =\mathcal{O}\left(h^{-2}\right)$.

Assuming that $y\left(x\right)$ has in $[0,x_N]$ bounded derivatives up to the necessary order, from (24) and the vector $L\left(h\right)$ of local truncation errors, we can obtain:

$$\begin{array}{cc}\hfill \parallel E\parallel \le & \parallel \left(M^{-1}\right)\parallel \parallel L\left(h\right)\parallel \hfill \\ \hfill =& \mathcal{O}\left(h^{-2}\right)\mathcal{O}\left(h^7\right)\hfill \\ \hfill \le & K{h}^5.\hfill \end{array}$$

We have shown that the global method exhibits a fifth-order convergence. Nevertheless, in view of the form of the vector $L\left(h\right)$, we see that, assuming the sufficient smoothness of the solution, at the mesh points we obtain a superconvergence order (see Ascher et al. [25]):

$$|e_j|=|y\left(x_j\right)-y_j|\le |\mathcal{O}\left(h^{-2}\right)\left|\right|\mathcal{O}\left(h^8\right)|\le K{h}^6,j=1,2,\dots N.$$

Therefore, the proposed method is convergent, providing sixth-order approximations. □

4. Implementation Issues

The PHNT is implemented in a block unification mode. We rewrite the systems in (20) as $\mathit{F}\left(\mathit{y}\right)=\mathbf{0}$ and the unknowns as:

$$\tilde{\mathbf{U}}=\bigcup {\left\{y_j\right\}}_{j=1,\dots ,N-1}\bigcup {\left\{y_j^{\prime }\right\}}_{j=0,\dots ,N}\bigcup {\left\{y_{j+u},y_{j+v},y_{j+w},y_{j+u}^{\prime },y_{j+v}^{\prime },y_{j+w}^{\prime }\right\}}_{j=0,1,\dots ,N-1}.$$

Then, we use Modified Newton’s method (MNM) to solve non-linear equations, since the PHNT is an implicit scheme. The MNM is given by:

$${\tilde{\mathbf{U}}}^{\mathbf{i}+\mathbf{1}}={\tilde{\mathbf{U}}}^{\mathbf{i}}-{\left({\mathbf{J}}^{\mathbf{i}}\right)}^{-\mathbf{1}}{\mathbf{F}}^{\mathbf{i}},$$

where $\mathbf{J}$ represents the jacobian matrix of $\mathbf{F}$. The starting values for using MNM for solving the systems given in (12)–(14) for each iteration are taken as those provided by the linear interpolation obtained throughout the boundary values, while the stopping criterion considers a maximum number of 100 iterations and an error between two successive approximations of less than ${10}^{-16}$.

We enumerate and summarize how the PHNT is utilized to give numerical solutions to physical models and catalytic diffusion–reaction problems as follows:

• Let us take $N>0\in \mathbb{N}$, and define $h=\frac{x_N-x_0}{N}$ to generate the partition:

  $$P_N=\bigcup {\left\{x_j\right\}}_{j=0,1,\dots ,N}\bigcup {\left\{x_{j+k}\right\}}_{k=u,v,w;j=0,1,\dots ,N-1}.$$
• Using Equations (12)–(14) and Equations (9)–(11) for $n=1,\dots ,N-1$, we can form a system of equations with variables:

  $$\left\{y_u,y_v,y_w,y_u^{\prime },y_v^{\prime },y_w^{\prime }\right\}\bigcup {\left\{y_j\right\}}_{j=1,\dots ,N-1}\bigcup {\left\{y_j^{\prime }\right\}}_{j=0,\dots ,N}\bigcup {\left\{y_{j+u},y_{j+v},y_{j+w},y_{j+u}^{\prime },y_{j+v}^{\prime },y_{j+w}^{\prime }\right\}}_{j=1,\dots ,N-1}.$$
• We make just one block matrix equation by joining all the equations generated in the previous step of the partition $P_N$ with the given boundary conditions.
• We solve the single block matrix equation simultaneously to obtain the approximate solutions for the SBVP on the whole interval $\left[x_0,x_N\right]$.

5. Numerical Illustrations

This section presents the numerical outcomes and discussion of the proposed PHNT for the solution of the singular physical models and catalytic diffusion–reaction problems of the form $\left(1\right)$. The accuracy of the PHNT is measured by utilizing the following formulas:

$$ABER=\parallel y\left(x_j\right)-y_j\parallel ,MAXABER=\underset{j=0,1,\dots ,N}{max}\parallel y\left(x_j\right)-y_j\parallel ,$$

where ABER denotes the absolute error at the considered node, MAXABER is the maximum absolute error along the considered interval, $y\left(x_j\right)$ is the theoretical solution, and $y_j$ is the approximate solution provided by the PHNT.

5.1. Example 1

We firstly consider the following scalar Lane–Emden singular equation (SCLSE), which corresponds to the reaction–diffusion process in a spherical permeable catalyst as reported in [3,26],

$$\begin{array}{c}\hfill y^{″}\left(x\right)+\frac{2}{x}{y}^{\prime }\left(x\right)-{\varphi }^{2}y{\left(x\right)}^n=0,y\left(1\right)=1\left(\mathrm{at}\mathrm{the}\mathrm{catalyst}\mathrm{surface}\right),\\ \hfill y^{\prime }\left(0\right)=0\left(\mathrm{at}\mathrm{the}\mathrm{centre}\mathrm{of}\mathrm{the}\mathrm{catalyst}\mathrm{surface}\right),x\in [0,1].\end{array}$$

(25)

The general analytical solution of problem (25) is unknown, but its solution for $n=1,$ is given by $y\left(x\right)=\frac{sinh\left(x\varphi \right)}{xsinh\left(\varphi \right)}$, where $\varphi$ is the Thiele modulus and ${\varphi }^2=\frac{\mathrm{reaction}\mathrm{rate}\mathrm{at}\mathrm{the}\mathrm{catalyst}\mathrm{surface}}{\mathrm{diffusion}\mathrm{rate}\mathrm{at}\mathrm{the}\mathrm{catalyst}\mathrm{pores}}$.

Table 1. Comparison of $PHNT$ and the exact solution on test 1 with $h=\frac{1}{10},n=1,\varphi =5$.

x	$\mathit{PHNT}$	Exact Solution	ABER
$0.1$	$0.07022543735304641$	$0.07022543922779090$	$1.87474\times {10}^{-9}$
$0.2$	$0.07918802960817403$	$0.07918802869127974$	$9.16894\times {10}^{-10}$
$0.3$	$0.09565082522028227$	$0.09565082330512954$	$1.91515\times {10}^{-9}$
$0.4$	$0.12219351619982496$	$0.12219351358270766$	$2.61712\times {10}^{-9}$
$0.5$	$0.16307123522003290$	$0.16307123192997786$	$3.29006\times {10}^{-9}$
$0.6$	$0.22500992040856140$	$0.22500991644891966$	$3.95964\times {10}^{-9}$
$0.7$	$0.31848116606908733$	$0.31848116156439320$	$4.50469\times {10}^{-9}$
$0.8$	$0.45971591487821106$	$0.45971591027923203$	$4.59898\times {10}^{-9}$
$0.9$	$0.67387038375554700$	$0.67387038020431460$	$3.55123\times {10}^{-9}$
1	$1.00000000000000000$	$1.0000000000000000$	$0.00000$

presents the numerical results with the proposed method, showing that they are very close to the theoretical solution available for $n=1$. The CPU time with the PHNT for the value of $N=10$ in Table 1 is $0.2188$ s.

To analyze the impact of the Thiele modulus $\left(\varphi \right)$ on the concentration profile $\left(y\right(x\left)\right)$, we also considered other values of $\varphi$ and n. Figure 1 displays the numerical outcomes for various values of $\varphi$ and n. We observed that in Figure 1, the concentration profile increases when $\varphi$ diminishes.

5.2. Example 2

As a second test problem, we consider the non-homogeneous SCLSE, which corresponds to the physical model problem 2 in [4]:

$$\begin{array}{c}\hfill y^{″}\left(x\right)+\frac{1}{x}{y}^{\prime }\left(x\right)+y\left(x\right)=4-9x+x^2-x^3,y\left(0\right)=0,y\left(1\right)=0\end{array}$$

(26)

where the exact solution is given by $y\left(x\right)=x^2-x^3$.

We have applied PHNT to test problem 2; Figure 2 presents the numerical solution of PHNT, which is very close to the exact solution for problem 2. Figure 3 shows a graphical representation of the absolute errors (ABER) for different values of x. The CPU time with the PHNT for the value of $N=10$ in Figure 2 and Figure 3 is $0.0625$ s.

5.3. Example 3

Let us consider the non-linear SCLSE, which corresponds to the physical model problem of thermal explosion in cylindrical vessel reported by Thula and Roul [3] and Roul et al. [14]:

$$\begin{array}{c}\hfill y^{″}\left(x\right)+\frac{1}{x}{y}^{\prime }\left(x\right)=exp\left(y\left(x\right)\right),y\left(0\right)=0,y^{\prime }\left(1\right)=0,\end{array}$$

(27)

where the analytical solution is given by $y\left(x\right)=2log\left(\frac{2\sqrt{6}+1-5}{\left(2\sqrt{6}-5\right)x^2+1}\right)$.

Test problem 3 is numerically solved using the new $PHNT$ for different values of step size $\left(h\right)$. The numerical results and the comparisons of MAXABER and ABER between $PHNT$ and the methods in [3,14] are abridged in

Table 2. Maximum absolute errors (MAXABER) for test problem 3.

h	Method	$\mathit{MAXABER}$	$\mathit{ROC}$
$\frac{1}{8}$	$PHNT$	$1.91969\times {10}^{-10}$
$\frac{1}{8}$	Method in [3]	$1.01400\times {10}^{-8}$
$\frac{1}{8}$	Method in [14]	$8.53810\times {10}^{-10}$
$\frac{1}{16}$	$PHNT$	$2.99397\times {10}^{-12}$	$6.00$
$\frac{1}{16}$	Method in [3]	$5.80500\times {10}^{-10}$	$4.13$
$\frac{1}{16}$	Method in [14]	$2.19100\times {10}^{-11}$	$5.30$
$\frac{1}{32}$	$PHNT$	$4.77118\times {10}^{-14}$	$5.97$
$\frac{1}{32}$	Method in [3]	$3.49500\times {10}^{-11}$	$4.05$
$\frac{1}{32}$	Method in [14]	$3.92400\times {10}^{-13}$	$5.94$

.

We note that the $PHNT$ method presents a better performance compared with the techniques in [3,14]. Additionally, the CPU time in our proposed $PHNT$ to obtain the approximate solutions for problem 3 with step size $h=\frac{1}{16}$ is $0.4531\text{s}.$ In Table 2, we have included the numerical rate of convergence (ROC) with the following formula:

$$ROC\simeq -{log}_2\left(\frac{MAXABE{R}_h}{MAXABE{R}_{2h}}\right).$$

5.4. Example 4

As a test problem 4, we consider a non-linear homogeneous SCLSE, which corresponds to the physical model problem arising in chemistry and chemical kinetics. The formulation of heat and mass transfer within porous catalyst particles is reported by Ravikanth [27]:

$$\begin{array}{c}\hfill y^{″}\left(x\right)+\frac{\lambda }{x}{y}^{\prime }\left(x\right)-{\varphi }^{2}y\left(x\right)exp\left(\frac{rs(1-y\left(x\right))}{c(1-y\left(x\right))+1}\right)=0,y\left(0\right)=0,y^{\prime }\left(1\right)=0,\end{array}$$

(28)

where the analytical solution is unknown.

We solved problem 4 with the PHNT scheme for $\varphi =r=s=c=1$. The PHNT approximate solutions for $\lambda =2$ and $\lambda =4$ for different values of x are plotted in Figure 4.

Table 3. Comparison of approximation solutions for test problem 4 for $\lambda =2$.

x	$\mathit{PHNT},\mathit{h}=\frac{1}{10}$	Method in [13], $\mathit{h}=\frac{1}{10}$	Method in [27], $\mathit{h}=\frac{1}{50}$
$0.1$	$0.8383648968983356$	$0.838364878696000$	$0.83836491959750$
$0.2$	$0.8431842515107244$	$0.843184233589000$	$0.84318428772800$
$0.3$	$0.8512302074850839$	$0.851230190133000$	$0.85123026453741$
$0.4$	$0.8625224114697174$	$0.862522394979000$	$0.86252249405263$
$0.5$	$0.8770865272403503$	$0.877086511985000$	$0.87708663616863$
$0.6$	$0.8949523006678164$	$0.894952287104000$	$0.89495243141739$
$0.7$	$0.9161509382771176$	$0.916150926969000$	$0.91615107927542$
$0.8$	$0.9407117001197036$	$0.940711691749000$	$0.94071183074544$
$0.9$	$0.9686575885218047$	$0.968657583887000$	$0.96865767679753$
$1.0$	$1.0000000000000000$	$1.0000000000000000$	$1.000000000000000$

contains the numerical results provided by PHNT. We also used, for comparison, the cubic spline method (CSM) and the modified Adomian decomposition technique (MADT) in [13,27], with the PHNT being significantly better. The CPU time used by PHNT for $N=10$ with the specifications in Table 3 is $0.5938$ s.

6. Conclusions

This work presented a reliable PHNT approach for solving the scalar and system of SBVP of Lane-Emden type in various physical models and chemical kinetics. We employed a set of optimized hybrid block formulas given in (9)–(11) which are combined with an appropriate starting algorithm in (12)–(14) specifically designed to cope with the singularity at the beginning of the integration interval of the considered problem. Four real-world model problems in applied sciences and engineering are solved numerically to show the strength and quality of the proposed PHNT. The obtained approximate results in Table 1, Table 2 and Table 3 and Figure 1, Figure 2, Figure 3 and Figure 4 accentuate the effectiveness of the new methodology. An interesting question to be addressed in future investigations is how to select the optimal stepsizes. Notice that we have used the same step size for both groups of formulas, but we could have chosen one stepsize for the formulas in (9)–(11) and another for the formulas in (12)–(14). How to obtain the optimal values for these step sizes is an open question.