An Improved Symmetric Numerical Approach for Systems of Second-Order Two-Point BVPs

Abstract

This study deals with the numerical solution of a class of linear systems of second-order boundary value problems (BVPs) using a new symmetric cubic B-spline method (NCBM). This is a typical cubic B-spline collocation method powered by new approximations for second-order derivatives. The flexibility and high order precision of B-spline functions allow them to approximate the answers. These functions have a symmetrical property. The new second-order approximation plays an important role in producing more accurate results up to a fifth-order accuracy. To verify the proposed method’s accuracy, it is tested on three linear systems of ordinary differential equations with multiple step sizes. The numerical findings by the present method are quite similar to the exact solutions available in the literature. We discovered that when the step size decreased, the computational errors decreased, resulting in better precision. In addition, details of maximum errors are investigated. Moreover, simple implementation and straightforward computations are the main advantages of the offered method. This method yields improved results, even if it does not require using free parameters. Thus, it can be concluded that the offered scheme is reliable and efficient.

1. Introduction

Most problems arising from scientific and engineering applications, especially applications in geodesics, are boundary value problems (BVPs), which are much more difficult to solve than initial value problems (IVPs). Since it is generally difficult to find closed-form solutions for BVPs, many researchers have attempted to develop methods to find approximate or numerical solutions for BVPs. Well-known methods involve the shooting method [1], finite difference methods [2,3,4] and spectral methods [5,6,7]. In some real-life situations, the shooting method produces numerically sensitive systems of algebraic equations, which must be solved using other numerical methods [8].

In the present paper, we consider the following system of linear two-point second-order BVPs:

(1)$\begin{array}{cc}\hfill & w^{″}\left(x\right)+p_1\left(x\right)w^{\prime }\left(x\right)+p_2\left(x\right)w\left(x\right)+p_3\left(x\right)z^{″}\left(x\right)+p_4\left(x\right)z^{\prime }\left(x\right)+p_5\left(x\right)z\left(x\right)=f_1\left(x\right),\hfill \\ \hfill & z^{″}\left(x\right)+q_1\left(x\right)z^{\prime }\left(x\right)+q_2\left(x\right)z\left(x\right)+q_3\left(x\right)w^{″}\left(x\right)+q_4\left(x\right)w^{\prime }\left(x\right)+q_5\left(x\right)w\left(x\right)=f_2\left(x\right),\hfill \end{array}$

with boundary conditions

(2)$\begin{array}{c}\hfill w\left(0\right)=w\left(1\right)=z\left(0\right)=z\left(1\right)=0,\end{array}$

where $0\le x\le 1$. In particular, $f_1\left(x\right)$ and $f_2\left(x\right)$ are given functions, and $p_i$ and $q_i$ with $i=1,2,3,4,5$ are continuous and sufficiently smooth functions on the interval $I=[0,1]$. Theorems that systematically list the existence and uniqueness of the problem solutions of (1) and (2) have been studied in [9]. In recent times, applications of linear and non-linear systems of two-point boundary value problems can be found in economics, biology, physics and mathematics. For instance, Nikooeinejad et al. [10] obtained the approximate solution of two-point BVPS for four applications of differential games in economics and management science using a combined numerical algorithm. In biology, the Shortley–Weller scheme has been implemented for a two-point boundary value problem. This numerical scheme was later applied to investigate tumor growth problems in heterogeneous microenvironments [11]. On the other hand, the application of two-point boundary value problems has been addressed in the problem of calculating rocket trajectories in the atmosphere [12].

Several researchers have investigated the linear and non-linear systems of second-order boundary value problems and produced various efficient and accurate numerical methods. These methods include the Laplace homotopy analysis [13,14], continuous genetic algorithm method [15], sinc collocation method [16,17], He’s homotopy perturbation method [18], reproducing kernel space method [19], multistage optimal homotopy asymptotic method [20], variational iteration method (VIM) [21] and Chebyshev finite difference method [22].

Researchers have been interested in the families of B-splines for their potential to approximate the solution of BVPs accurately and efficiently. B-spline methods have several attractive features and flexibility that make them useful in numerical computation to solve BVPs [23]. For example, the B-spline is the smoothest interpolation function compared to other piecewise polynomial interpolation functions [24]. Moreover, B-splines have small local support properties. In recent years, the cubic B-spline collocation method captured the attention of some researchers to solve partial differential equations [25], fractional differential equations [26], fractional partial differential equations [27], etc.

This study focuses on finding the solutions of two-point BVPs using the cubic B-spline method. Bickley was the first to explore cubic splines to approximate the solutions of two-point BVPs [28]. Later, Albasiny and Hoskins enhanced Bickley’s work by solving the two-point BVPs using a tri-diagonal matrix of coefficients [29]. Since then, several researchers have earned more interest in employing spline functions for solving BVPs [30,31,32,33]. Caglar et al. in [34] evaluated the two-point BVPs solutions using the cubic B-spline basis function. Hamid et al. [35,36] considered the ECBM and cubic trigonometric B-spline method for the solution of linear two-point BVPs. Apart from that, Heilat and Ismail [37] used a hybrid cubic B-spline method to evaluate the solutions of non-linear two-point boundary value problems. Recently, a hybrid cubic B-spline method with an optimized parameter was used by Heilat et al. [38] to solve linear two-point BVPs.

The linear system of second-order BVPs has gained attention from Caglar and Caglar [39] and Heilat et al. [40]. They represented the cubic B-spline method (CBM) and ECBM, respectively. The ECBM involved two parameters in boosting the flexibility of the spline curve. Based on the investigation, the ECBM is the best compared to the CBM, He’s homotopy perturbation method [18], Laplace homotopy analysis method [13], reproducing kernel [19] and sinc-collocation method [17]. In recent years, Zhang and Niu [41] found the approximate solution of second-order BVPs using a Lobatto-reproducing kernel and declared the method has high precision accuracy in different spaces.

The new symmetric cubic B-spline method (NCBM) was first studied by Lang and Xu in [42] to solve non-linear second-order BVPs with two dependent variables. The NCBM is a typical CBM, equipped with a new second-order derivative approximation. Then, Iqbal et al. [43] explored the NCBM for solving several third-order Emden–Flower type equations. A year after that, Wasim et al. [44] extended the NCBM and proposed the new extended cubic B-spline method (NECBM) for solving the class of second-order singular BVPs. Moreover, the nonlinear third-order Korteweg–de Vries equations were solved by Abbas et al. [45] using the NCBM to approximate the solutions. Later, Nazir et al. [46] improved the method to a new quintic B-spline approximation technique as a method to approximate the numerical solution of the Boussinesq equation. Recently, Nazir et al. [47] implemented the NCBM for the numerical solutions of coupled viscous Burgers equations.

Thus, motivated by all these works, we aim to figure out whether the proposed method, the NCBM, can perform much better in solving the linear system of two-point second-order BVPs. The rest of this paper is as follows. In Section 2, the typical definition of cubic B-spline basis functions is described. Then, Section 3 presents the descriptions of the numerical method. The convergence analysis of the method has been proven in Section 4. The numerical results and their discussion are summarized in Section 5. Finally, the paper ends in Section 6 with a brief conclusion.

2. Preliminary Concepts

This section describes the classical cubic B-spline approximation and the new second-order approximation invented by Lang and Xu [42]. Let the finite interval $[a,b]$, where $a=x_0<\dots <x_N=b$ is divided into uniform partitions with a mesh point $x_i=x_0+ih,i=0\left(1\right)N$ using a step size $h=\frac{b-a}{N},N\in {\mathbb{Z}}^{+}$. The typical cubic B-spline basis function is defined as [34].

(3)$B_i\left(x\right)=\left\{\begin{array}{cc}{(x-x_{i-2})}^3,\hfill & x\in [x_{i-2},x_{i-1}],\hfill \\ h^3+3h^2(x-x_{i-1})+3h{(x-x_{i-1})}^2-3{(x-x_{i-1})}^3,\hfill & x\in [x_{i-1},x_i],\hfill \\ h^3+3h^2(x_{i+1}-x)+3h{(x_{i+1}-x)}^2-3{(x_{i+1}-x)}^3,\hfill & x\in [x_i,x_{i+1}],\hfill \\ {(x_{i+2}-x)}^3,\hfill & x\in [x_{i+1},x_{i+2}],\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$

where $i=-1\left(1\right)N+1$. The cubic B-spline holds the geometric invariability, convex hull property and symmetry [48]. For a sufficiently smooth function $w\left(x\right)$ and $z\left(x\right)$, there exist a unique third-degree spline $W\left(x\right)$ and $Z\left(x\right)$ that satisfies the prescribed interpolating conditions given by

$\begin{array}{cc}\hfill & W\left(x_i\right)=w\left(x_i\right),i=0,1,2,...,N,\hfill \\ \hfill & W^{\prime }\left(a\right)=w^{\prime }\left(a\right),W^{\prime }\left(b\right)=w^{\prime }\left(b\right),\hfill \\ \hfill & W^{″}\left(a\right)=w^{″}\left(a\right),W^{″}\left(b\right)=w^{″}\left(b\right),\hfill \end{array}$

and

$\begin{array}{cc}\hfill & Z\left(x_i\right)=z\left(x_i\right),i=0,1,2,...,N,\hfill \\ \hfill & Z^{\prime }\left(a\right)=z^{\prime }\left(a\right),Z^{\prime }\left(b\right)=z^{\prime }\left(b\right),\hfill \\ \hfill & Z^{″}\left(a\right)=z^{″}\left(a\right),Z^{″}\left(b\right)=z^{″}\left(b\right),\hfill \end{array}$

in which

(4)$W\left(x\right)=\sum _{i=-1}^{N+1}{\sigma }_i{B}_i\left(x\right),$

(5)$Z\left(x\right)=\sum _{i=-1}^{N+1}{\eta }_i{B}_i\left(x\right),$

where ${\sigma }_i$ and ${\eta }_i$ are unknown real coefficients to be computed. The values of $B_i\left(x\right)$ and the first and second derivatives $B_i^{\prime }\left(x\right)$ and $B_i^{″}\left(x\right)$ at mesh point $x_i$ are tabulated in Table 1. From (4), (5) and Table 1, the cubic B-spline approximations $W\left(x_j\right)$, $W^{\prime }\left(x_j\right)$, $Z\left(x_j\right)$ and $Z^{\prime }\left(x_j\right)$ can be simplified as follows:

(6)$\begin{array}{cc}\hfill W_j& =\sum _{i=j-1}^{j+1}{\sigma }_i{B}_i\left(x\right)=\frac{1}{6}({\sigma }_{j-i}+4{\sigma }_j+{\sigma }_{j+1}),\hfill \end{array}$

(7)$\begin{array}{cc}\hfill s_{j\left(x\right)}& =\sum _{i=j-1}^{j+1}{\sigma }_i{B}_i^{\prime }\left(x\right)=\frac{1}{2h}(-{\sigma }_{j-i}+{\sigma }_{j+1}),\hfill \end{array}$

(8)$\begin{array}{cc}\hfill Z_{j\left(x\right)}& =\sum _{i=j-1}^{j+1}{\eta }_i{B}_i\left(x\right)=\frac{1}{6}({\eta }_{j-i}+4{\eta }_j+{\eta }_{j+1}),\hfill \end{array}$

(9)$\begin{array}{cc}\hfill r_{j\left(x\right)}& =\sum _{i=j-1}^{j+1}{\eta }_i{B}_i^{\prime }\left(x\right)=\frac{1}{2h}(-{\eta }_{j-i}+{\eta }_{j+1}).\hfill \end{array}$

The second derivatives, $W^{″}\left(x_j\right)$ and $Z^{″}\left(x_j\right)$ can be simplified as $S_j$ and $R_j$, respectively. Subsequently, the new approximation for second-order derivatives can be represented as follows [42,49]:

(10)$S_j=\frac{1}{12h^2}\left\{\begin{array}{cc}14{\sigma }_{j-1}-33{\sigma }_j+28{\sigma }_{j+1}-14{\sigma }_{j+2}+6{\sigma }_{j+3}-{\sigma }_{j+4},\hfill & \mathrm{for}j=0,\hfill \\ {\sigma }_{j-2}+8{\sigma }_{j-1}-18{\sigma }_j+8{\sigma }_{j+1}+{\sigma }_{j+2},\hfill & \mathrm{for}j=1,\dots ,N-1,\hfill \\ -{\sigma }_{n-4}+6{\sigma }_{n-3}-14{\sigma }_{n-2}+28{\sigma }_{n-1}-33{\sigma }_n+14{\sigma }_{n+1},\hfill & \mathrm{for}j=N,\hfill \end{array}\right.$

and

(11)$R_j=\frac{1}{12h^2}\left\{\begin{array}{cc}14{\eta }_{j-1}-33{\eta }_j+28{\eta }_{j+1}-14{\eta }_{j+2}+6{\eta }_{j+3}-{\eta }_{j+4},\hfill & \mathrm{for}j=0,\hfill \\ {\eta }_{j-2}+8{\eta }_{j-1}-18{\eta }_j+8{\eta }_{j+1}+{\eta }_{j+2},\hfill & \mathrm{for}j=1,\dots ,N-1,\hfill \\ -{\eta }_{n-4}+6{\eta }_{n-3}-14{\eta }_{n-2}+28{\eta }_{n-1}-33{\eta }_n+14{\eta }_{n+1},\hfill & \mathrm{for}j=N.\hfill \end{array}\right.$

We note that this NCBM has up to a fifth-order accuracy [49].

3. Implementation of the Method

In this section, we extended Caglar’s work [39] by solving the linear system two-point second-order BVP and adopted the new second-order approximation.

Discretizing (1) at the knot $x_j$ gives the following expression:

(12)$\begin{array}{cc}\hfill & S_{k+1}\left(x_j\right)+p_1\left(x_j\right)s_{k+1}\left(x_j\right)+p_2\left(x_j\right)W_{k+1}\left(x_j\right)+p_3\left(x_j\right)R_{k+1}\left(x_j\right)+p_4\left(x_j\right)r_{k+1}\left(x_j\right)\hfill \\ \hfill & +p_5\left(x_j\right)Z_{k+1}\left(x_j\right)=f_{1k}\left(x_j\right),\hfill \end{array}$

(13)$\begin{array}{cc}\hfill & R_{k+1}\left(x_j\right)+q_1\left(x_j\right)r_{k+1}\left(x_j\right)+q_2\left(x_j\right)Z_{k+1}\left(x_j\right)+q_3\left(x_j\right)S_{k+1}\left(x_j\right)+q_4\left(x_j\right)s_{k+1}\left(x_j\right)\hfill \\ \hfill & +q_5\left(x_j\right)W_{k+1}\left(x_j\right)=f_{2k}\left(x_j\right),\hfill \end{array}$

where $k=0,1,2,\dots$. By substituting (6)–(11) into (12) for $j=0,1,2,\dots$, $N-1,N$, we obtain the following equation:

For $j=0$,

(14)$\begin{array}{cc}\hfill & \left(\frac{14{\sigma }_{-1}-33{\sigma }_0+28{\sigma }_1-14{\sigma }_2+6{\sigma }_3-{\sigma }_4}{12h^2}\right)\hfill \\ \hfill & +p_1\left(x_0\right)\left(\frac{-{\sigma }_{-1}+{\sigma }_1}{2h}\right)+p_2\left(x_0\right)\left(\frac{{\sigma }_{-1}+4{\sigma }_0+{\sigma }_1}{6}\right)\hfill \\ \hfill & +p_3\left(x_0\right)\left(\frac{14{\eta }_{-1}-33{\eta }_0+28{\eta }_1-14{\eta }_2+6{\eta }_3-{\eta }_4}{12h^2}\right)\hfill \\ \hfill & +p_4\left(x_0\right)\left(\frac{-{\eta }_{-1}+{\eta }_1}{2h}\right)+p_5\left(x_0\right)\left(\frac{{\eta }_{-1}+4{\eta }_0+{\eta }_1}{6}\right)=f_{k1}\left(x_0\right).\hfill \end{array}$

For $j=1,2,\dots ,N-1$,

(15)$\begin{array}{cc}\hfill & \left(\frac{{\sigma }_{j-2}+8{\sigma }_{j-1}-18{\sigma }_j+8{\sigma }_{j+1}+{\sigma }_{j+2}}{12h^2}\right)\hfill \\ \hfill & +p_1\left(x_j\right)\left(\frac{-{\sigma }_{j-1}+{\sigma }_{j+1}}{2h}\right)+p_2\left(x_j\right)\left(\frac{{\sigma }_{j-1}+4{\sigma }_j+{\sigma }_{j+1}}{6}\right)\hfill \\ \hfill & +p_3\left(x_j\right)\left(\frac{{\eta }_{j-2}+8{\eta }_{j-1}-18{\eta }_j+8{\eta }_{j+1}+{\eta }_{j+2}}{12h^2}\right)\hfill \\ \hfill & +p_4\left(x_j\right)\left(\frac{-{\eta }_{j-1}+{\eta }_{j+1}}{2h}\right)+p_5\left(x_j\right)\left(\frac{{\eta }_{j-1}+4{\eta }_j+{\eta }_{j+1}}{6}\right)=f_{k1}\left(x_j\right).\hfill \end{array}$

For $j=N$,

(16)$\begin{array}{cc}\hfill & \left(\frac{-{\sigma }_{N-4}+6{\sigma }_{N-3}-14{\sigma }_{N-2}+28{\sigma }_{N-1}-33{\sigma }_N+14{\sigma }_{N+1}}{12h^2}\right)\hfill \\ \hfill & +p_1\left(x_N\right)\left(\frac{-{\sigma }_{N-1}+{\sigma }_{N+1}}{2h}\right)+p_2\left(x_N\right)\left(\frac{{\sigma }_{N-1}+4{\sigma }_N+{\sigma }_{N+1}}{6}\right)\hfill \\ \hfill & +p_4\left(x_N\right)\left(\frac{-{\eta }_{N-4}+6{\eta }_{N-3}-14{\eta }_{N-2}+28{\eta }_{N-1}-33{\eta }_N+14{\eta }_{N+1}}{12h^2}\right)\hfill \\ \hfill & +p_4\left(x_N\right)\left(\frac{-{\eta }_{N-1}+{\eta }_{N+1}}{2h}\right)+p_5\left(x_N\right)\left(\frac{{\eta }_{N-1}+4{\eta }_N+{\eta }_{N+1}}{6}\right)=f_{k1}\left(x_N\right).\hfill \end{array}$

By substituting (6)–(11) into (13) for $j=0,1,2,\dots$, $N-1,N$, we obtain the following equation:

For $j=0$,

(17)$\begin{array}{cc}\hfill & \left(\frac{14{\eta }_{-1}-33{\eta }_0+28{\eta }_1-14{\eta }_2+6{\eta }_3-{\eta }_4}{12h^2}\right)\hfill \\ \hfill & +q_1\left(x_0\right)\left(\frac{-{\eta }_{-1}+{\eta }_1}{2h}\right)+q_2\left(x_0\right)\left(\frac{{\eta }_{-1}+4{\eta }_0+{\eta }_1}{6}\right)\hfill \\ \hfill & +q_3\left(x_0\right)\left(\frac{14{\sigma }_{-1}-33{\sigma }_0+28{\sigma }_1-14{\sigma }_2+6{\sigma }_3-{\sigma }_4}{12h^2}\right)\hfill \\ \hfill & +q_4\left(x_0\right)\left(\frac{-{\sigma }_{-1}+{\sigma }_1}{2h}\right)+q_5\left(x_0\right)\left(\frac{{\sigma }_{-1}+4{\sigma }_0+{\sigma }_1}{6}\right)=f_{k1}\left(x_0\right).\hfill \end{array}$

For $j=1,2,\dots ,N-1$,

(18)$\begin{array}{cc}\hfill & \left(\frac{{\eta }_{j-2}+8{\eta }_{j-1}-18{\eta }_j+8{\eta }_{j+1}+{\eta }_{j+2}}{12h^2}\right)\hfill \\ \hfill & +q_1\left(x_j\right)\left(\frac{-{\eta }_{j-1}+{\eta }_{j+1}}{2h}\right)+q_2\left(x_j\right)\left(\frac{{\eta }_{j-1}+4{\eta }_j+{\eta }_{j+1}}{6}\right)\hfill \\ \hfill & +q_3\left(x_j\right)\left(\frac{{\sigma }_{j-2}+8{\sigma }_{j-1}-18{\sigma }_j+8{\sigma }_{j+1}+{\sigma }_{j+2}}{12h^2}\right)\hfill \\ \hfill & +q_4\left(x_j\right)\left(\frac{-{\sigma }_{j-1}+{\sigma }_{j+1}}{2h}\right)+q_5\left(x_j\right)\left(\frac{{\sigma }_{j-1}+4{\sigma }_j+{\sigma }_{j+1}}{6}\right)=f_{k1}\left(x_j\right).\hfill \end{array}$

For $j=N$,

(19)$\begin{array}{cc}\hfill & \left(\frac{-{\eta }_{N-4}+6{\eta }_{N-3}-14{\eta }_{N-2}+28{\eta }_{N-1}-33{\eta }_N+14{\eta }_{N+1}}{12h^2}\right)\hfill \\ \hfill & +q_1\left(x_N\right)\left(\frac{-{\eta }_{N-1}+{\eta }_{N+1}}{2h}\right)+q_2\left(x_N\right)\left(\frac{{\eta }_{N-1}+4{\eta }_N+{\eta }_{N+1}}{6}\right)\hfill \\ \hfill & +q_3\left(x_N\right)\left(\frac{-{\sigma }_{N-4}+6{\sigma }_{N-3}-14{\sigma }_{N-2}+28{\sigma }_{N-1}-33{\sigma }_N+14{\sigma }_{N+1}}{12h^2}\right)\hfill \\ \hfill & +q_4\left(x_N\right)\left(\frac{-{\sigma }_{N-1}+{\sigma }_{N+1}}{2h}\right)+q_5\left(x_N\right)\left(\frac{{\sigma }_{N-1}+4{\sigma }_N+{\sigma }_{N+1}}{6}\right)=f_{k1}\left(x_N\right).\hfill \end{array}$

Consequently, we have $2N+2$ linear equations involving $2N+6$ unknowns. Thus, we need four additional equations, which can be obtained from the boundary conditions in (2) below:

$\begin{array}{cc}\hfill & {\sigma }_{-1}+4{\sigma }_0+{\sigma }_1=0,\hfill \\ \hfill & {\sigma }_{N-1}+4{\sigma }_N+{\sigma }_{N+1}=0,\hfill \\ \hfill & {\eta }_{-1}+4{\eta }_0+{\eta }_1=0,\hfill \\ \hfill & {\eta }_{N-1}+4{\eta }_N+{\eta }_{N+1}=0.\hfill \end{array}$

Hence, the above system will have the $(2N+6)\times (2N+6)$ dimensional matrix form that can be expressed as:

$AX=B,$

where matrix A is given by:

$A=\left[\begin{array}{ccc}{A}_1& |& A_2\\ \cdots & \cdots & \cdots \\ A_4& |& A_3\end{array}\right],$

$\begin{array}{ccc}\hfill X& =& {[{\sigma }_{-1},{\sigma }_0,{\sigma }_1,\dots ,{\sigma }_{N-1},{\sigma }_N,{\sigma }_{N+1},{\eta }_{-1},{\eta }_0,{\eta }_1,\dots ,{\eta }_{N-1},{\eta }_N,{\eta }_{N+1}]}^T\hfill \end{array}$

and

$\begin{array}{ccc}\hfill B& =& [0,12h^2{f}_1\left(x_0\right),12h^2{f}_1\left(x_1\right),\dots ,12h^2{f}_1\left(x_N\right),0,0,12h^2{f}_2\left(x_0\right),12h^2{f}_2\left(x_1\right),\hfill \\ & & \dots ,12h^2{f}_2\left(x_N\right),0{\left)\right]}^T.\hfill \end{array}$

The four sub-matrices $A_1$, $A_2$, $A_3$ and $A_4$ are represented as follows:

$A_1=\left(\begin{array}{c}\begin{array}{ccccccccccccc}1& 4& 1& 0& 0& 0& \cdots & 0& 0& 0& 0& 0& 0\\ {ϵ}_1& {ϵ}_2& {ϵ}_3& {ϵ}_4& {ϵ}_5& {ϵ}_6& \cdots & 0& 0& 0& 0& 0& 0\\ {\alpha }_{1,1}& {\alpha }_{2,1}& {\alpha }_{3,1}& {\alpha }_{4,1}& {\alpha }_{5,1}& 0& \cdots & 0& 0& 0& 0& 0& 0\\ 0& {\alpha }_{1,2}& {\alpha }_{2,2}& {\alpha }_{3,2}& {\alpha }_{4,2}& {\alpha }_{5,2}& \cdots & 0& 0& 0& 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 0& 0& \cdots & {\alpha }_{1,N-2}& {\alpha }_{2,N-2}& {\alpha }_{3,N-2}& {\alpha }_{4,N-2}& {\alpha }_{5,N-2}& 0\\ 0& 0& 0& 0& 0& 0& \cdots & 0& {\alpha }_{1,N-1}& {\alpha }_{2,N-1}& {\alpha }_{3,N-1}& {\alpha }_{4,N-1}& {\alpha }_{5,N-1}\\ 0& 0& 0& 0& 0& 0& \cdots & {\varphi }_1& {\varphi }_2& {\varphi }_3& {\varphi }_4& {\varphi }_5& {\varphi }_6\\ 0& 0& 0& 0& 0& 0& \cdots & 0& 0& 0& 1& 4& 1\end{array}\end{array}\right),$

$A_2=\left(\begin{array}{c}\begin{array}{ccccccccccccc}0& 0& 0& 0& 0& 0& \cdots & 0& 0& 0& 0& 0& 0\\ {ϵ}_7& {ϵ}_8& {ϵ}_9& {ϵ}_{10}& {ϵ}_{11}& {ϵ}_{12}& \cdots & 0& 0& 0& 0& 0& 0\\ {\alpha }_{6,1}& {\alpha }_{7,1}& {\alpha }_{8,1}& {\alpha }_{9,1}& {\alpha }_{10,1}& 0& \cdots & 0& 0& 0& 0& 0& 0\\ 0& {\alpha }_{6,2}& {\alpha }_{7,2}& {\alpha }_{8,2}& {\alpha }_{9,2}& {\alpha }_{10,2}& \cdots & 0& 0& 0& 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 0& 0& \cdots & {\alpha }_{6,N-2}& {\alpha }_{7,N-2}& {\alpha }_{8,N-2}& {\alpha }_{9,N-2}& {\alpha }_{10,N-2}& 0\\ 0& 0& 0& 0& 0& 0& \cdots & 0& {\alpha }_{6,N-1}& {\alpha }_{7,N-1}& {\alpha }_{8,N-1}& {\alpha }_{9,N-1}& {\alpha }_{10,N-1}\\ 0& 0& 0& 0& 0& 0& \cdots & {\varphi }_7& {\varphi }_8& {\varphi }_9& {\varphi }_{10}& {\varphi }_{11}& {\varphi }_{12}\\ 0& 0& 0& 0& 0& 0& \cdots & 0& 0& 0& 0& 0& 0\end{array}\end{array}\right),$

$A_3=\left(\begin{array}{c}\begin{array}{ccccccccccccc}1& 4& 1& 0& 0& 0& \cdots & 0& 0& 0& 0& 0& 0\\ {\beta }_1& {\beta }_2& {\beta }_3& {\beta }_4& {\beta }_5& {\beta }_6& \cdots & 0& 0& 0& 0& 0& 0\\ {\gamma }_{1,1}& {\gamma }_{2,1}& {\gamma }_{3,1}& {\gamma }_{4,1}& {\gamma }_{5,1}& 0& \cdots & 0& 0& 0& 0& 0& 0\\ 0& {\gamma }_{1,2}& {\gamma }_{2,2}& {\gamma }_{3,2}& {\gamma }_{4,2}& {\gamma }_{5,2}& \cdots & 0& 0& 0& 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 0& 0& \cdots & {\gamma }_{1,N-2}& {\gamma }_{2,N-2}& {\gamma }_{3,N-2}& {\gamma }_{4,N-2}& {\gamma }_{5,N-2}& 0\\ 0& 0& 0& 0& 0& 0& \cdots & 0& {\gamma }_{1,N-1}& {\gamma }_{2,N-1}& {\gamma }_{3,N-1}& {\gamma }_{4,N-1}& {\gamma }_{5,N-1}\\ 0& 0& 0& 0& 0& 0& \cdots & {\delta }_1& {\delta }_2& {\delta }_3& {\delta }_4& {\delta }_5& {\delta }_6\\ 0& 0& 0& 0& 0& 0& \cdots & 0& 0& 0& 1& 4& 1\end{array}\end{array}\right),$

and

$A_4=\left(\begin{array}{c}\begin{array}{ccccccccccccc}0& 0& 0& 0& 0& 0& \cdots & 0& 0& 0& 0& 0& 0\\ {\beta }_7& {\beta }_8& {\beta }_9& {\beta }_{10}& {\beta }_{11}& {\beta }_{12}& \cdots & 0& 0& 0& 0& 0& 0\\ {\gamma }_{6,1}& {\gamma }_{7,1}& {\gamma }_{8,1}& {\gamma }_{9,1}& {\gamma }_{10,1}& 0& \cdots & 0& 0& 0& 0& 0& 0\\ 0& {\gamma }_{6,2}& {\gamma }_{7,2}& {\gamma }_{8,2}& {\gamma }_{9,2}& {\gamma }_{10,2}& \cdots & 0& 0& 0& 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 0& 0& \cdots & {\gamma }_{6,N-2}& {\gamma }_{7,N-2}& {\gamma }_{8,N-2}& {\gamma }_{9,N-2}& {\gamma }_{10,N-2}& 0\\ 0& 0& 0& 0& 0& 0& \cdots & 0& {\gamma }_{6,N-1}& {\gamma }_{7,N-1}& {\gamma }_{8,N-1}& {\gamma }_{9,N-1}& {\gamma }_{10,N-1}\\ 0& 0& 0& 0& 0& 0& \cdots & {\delta }_7& {\delta }_8& {\delta }_9& {\delta }_{10}& {\delta }_{11}& {\delta }_{12}\\ 0& 0& 0& 0& 0& 0& \cdots & 0& 0& 0& 0& 0& 0\end{array}\end{array}\right).$

For $j=0,$

$\begin{array}{cccc}\hfill {ϵ}_1& =14-6p_1\left(x_0\right)h+2p_2\left(x_0\right)h^2,\hfill & \hfill {ϵ}_7& =14p_3\left(x_0\right)-6p_4\left(x_0\right)+2p_5\left(x_0\right)h^2,\hfill \\ \hfill {ϵ}_2& =-33+8p_2\left(x_0\right)h^2,\hfill & \hfill {ϵ}_8& =-33p_3\left(x_0\right)+8p_5\left(x_0\right)h^2,\hfill \\ \hfill {ϵ}_3& =28+6p_1\left(x_0\right)h+2p_2\left(x_0\right)h^2,\hfill & \hfill {ϵ}_9& =28p_3\left(x_0\right)+6p_4\left(x_0\right)h+2p_5\left(x_0\right)h^2,\hfill \\ \hfill {ϵ}_4& =-14,\hfill & \hfill {ϵ}_{10}& =-14p_3\left(x_0\right),\hfill \\ \hfill {ϵ}_5& =6,\hfill & \hfill {ϵ}_{11}& =6p_3\left(x_0\right),\hfill \\ \hfill {ϵ}_6& =-1,\hfill & \hfill {ϵ}_{12}& =-p_3\left(x_0\right),\hfill \end{array}$

$\begin{array}{cccc}\hfill {\beta }_1& =14-6q_1\left(x_0\right)h+2q_2\left(x_0\right)h^2,\hfill & \hfill {\beta }_7& =14q_3\left(x_0\right)-6q_4\left(x_0\right)h+2q_5\left(x_0\right)h^2,\hfill \\ \hfill {\beta }_2& =-33+8q_2\left(x_0\right)h^2;\hfill & \hfill {\beta }_8& =-33q_3\left(x_0\right)+8q_5\left(x_0\right)h^2,\hfill \\ \hfill {\beta }_3& =28+6q_1\left(x_0\right)h+2q_2\left(x_0\right)h^2,\hfill & \hfill {\beta }_9& =28q_3\left(x_0\right)+6q_4\left(x_0\right)h+2q_5\left(x_0\right)h^2,\hfill \\ \hfill {\beta }_4& =-14,\hfill & \hfill {\beta }_{10}& =-14q_3\left(x_0\right),\hfill \\ \hfill {\beta }_5& =6,\hfill & \hfill {\beta }_{11}& =6q_3\left(x_0\right),\hfill \\ \hfill {\beta }_6& =-1,\hfill & \hfill {\beta }_{12}& =-q_3\left(x_0\right).\hfill \end{array}$

For $j=1,2,\dots ,N-1,$

$\begin{array}{cccc}\hfill {\alpha }_{1,j}& =1,\hfill & \hfill {\alpha }_{6,j}& =p_3\left(x_j\right),\hfill \\ \hfill {\alpha }_{2,j}& =8-6p_1\left(x_j\right)h+2p_2\left(x_j\right)h^2,\hfill & \hfill {\alpha }_{7,j}& =8p_3\left(x_j\right)-6p_4\left(x_j\right)h+2p_5\left(x_j\right)h^2,\hfill \\ \hfill {\alpha }_{3,j}& =-18+8p_2\left(x_j\right)h^2,\hfill & \hfill {\alpha }_{8,j}& =-18p_3\left(x_j\right)+8p_5\left(x_j\right)h^2,\hfill \\ \hfill {\alpha }_{4,j}& =8+6p_1\left(x_j\right)h+2p_2\left(x_j\right)h^2,\hfill & \hfill {\alpha }_{9,j}& =8p_3\left(x_j\right)+6p_4\left(x_j\right)h+2p_5\left(x_j\right)h^2,\hfill \\ \hfill {\alpha }_{5,j}& =1,\hfill & \hfill {\alpha }_{10,j}& =p_3\left(x_j\right),\hfill \end{array}$

$\begin{array}{cccc}\hfill {\gamma }_{1,j}& =1,\hfill & \hfill {\gamma }_{6,j}& =q_3\left(x_j\right),\hfill \\ \hfill {\gamma }_{2,j}& =8-6q_1\left(x_j\right)h+2q_2{h}^2,\hfill & \hfill {\gamma }_{7,j}& =8q_3\left(x_j\right)-6q_4\left(x_j\right)h+2q_5\left(x_j\right)h^2,\hfill \\ \hfill {\gamma }_{3,j}& =-18+8q_2\left(x_j\right)h^2,\hfill & \hfill {\gamma }_{8,j}& =-18q_3\left(x_j\right)+8q_5\left(x_j\right)h^2,\hfill \\ \hfill {\gamma }_{4,j}& =8+6q_{(}{x}_j)h+2q_2\left(x_j\right)h^2,\hfill & \hfill {\gamma }_{9,j}& =8q_3\left(x_j\right)+6q_4\left(x_j\right)h+2q_5\left(x_j\right)h^2,\hfill \\ \hfill {\gamma }_{5,j}& =1,\hfill & \hfill {\gamma }_{10,j}& =q3\left(x_j\right).\hfill \end{array}$

For $j=N$,

$\begin{array}{cccc}\hfill {\varphi }_1& =-1,\hfill & \hfill {\varphi }_7& =-p_3\left(x_N\right),\hfill \\ \hfill {\varphi }_2& =6,\hfill & \hfill {\varphi }_8& =6p_3\left(x_N\right),\hfill \\ \hfill {\varphi }_3& =-14,\hfill & \hfill {\varphi }_9& =-14p_3\left(x_N\right),\hfill \\ \hfill {\varphi }_4& =28-6p_1\left(x_N\right)h+2p_2\left(x_N\right)h^2,\hfill & \hfill {\varphi }_{10}& =28p_3\left(x_N\right)-6p_4\left(x_N\right)h+2p_5\left(x_N\right)h^2,\hfill \\ \hfill {\varphi }_5& =-33+8p_2\left(x_N\right)h^2,\hfill & \hfill {\varphi }_{11}& =-33p_3\left(x_N\right)+8p_5\left(x_N\right)h^2,\hfill \\ \hfill {\varphi }_6& =14+6p_1\left(x_N\right)h+2p_2\left(x_N\right)h^2,\hfill & \hfill {\varphi }_{12}& =14p_3\left(x_N\right)h+6p_4\left(x_N\right)h+2p_5\left(x_N\right)h^2,\hfill \end{array}$

$\begin{array}{cccc}\hfill {\delta }_1& =-1,\hfill & \hfill {\delta }_7& =-q_3\left(x_N\right),\hfill \\ \hfill {\delta }_2& =6,\hfill & \hfill {\delta }_8& =6q_3\left(x_N\right),\hfill \\ \hfill {\delta }_3& =-14,\hfill & \hfill {\delta }_9& =-14q_3\left(x_N\right),\hfill \\ \hfill {\delta }_4& =28-6q_1\left(x_N\right)h+2q_2\left(x_N\right)h^2,\hfill & \hfill {\delta }_{10}& =28q_3\left(x_N\right)-6q_4\left(x_N\right)h+2q_5\left(x_N\right)h^2,\hfill \\ \hfill {\delta }_5& =-33+8q_2\left(x_N\right)h^2,\hfill & \hfill {\delta }_{11}& =-33q_3\left(x_N\right)+8q_5\left(x_N\right)h^2,\hfill \\ \hfill {\delta }_6& =14+6q_1\left(x_N\right)h+2q_2\left(x_N\right)h^2,\hfill & \hfill {\delta }_{12}& =14q_3\left(x_N\right)+6q_4\left(x_N\right)h+2q_5\left(x_N\right)h^2.\hfill \end{array}$

Since the matrix A is a banded matrix, the system of linear equations is solved using a generalization of the Thomas algorithm. This method has been proposed in [50]. Matlab R2018a running on an Intel(R) CORE(TM) i7-1165G7 CPU 1.30 GHz processor, 8.00 GB RAM, was used to execute the numerical computations.

4. Convergence Analysis

In this section, we will prove the order of convergence of our method.

Theorem 1.

Let $p_i\left(x\right)\in C^2[0,1]$, where $i=1,2,3,4,5$ are continuous and sufficiently smooth functions. Then, let $\tilde{w}$ be the known exact solution of the boundary value problems (1), (2) and also $\tilde{m}$ be the cubic B-spline approximation to $\tilde{w}$. Thus, the uniform error is stated by

${\parallel \tilde{w}-\tilde{m}\parallel }_{\infty }\le \varphi h^2.$

Proof. 

Let $\tilde{w}$ be the exact solution of the boundary value problems (1), (2) and $\tilde{m}$ be the cubic B-spline approximation to $\tilde{w}$ given by:

(20)$\tilde{w}=\tilde{m}=\sum _{i=-1}^{N+1}\widetilde{{\sigma }_i}{B}_i{\left(x\right)}_i$

where

(21)$\tilde{\sigma }=\widetilde{{\sigma }_i}={[{\tilde{\sigma }}_{-1},{\tilde{\sigma }}_0,\dots ,{\tilde{\sigma }}_{N+1}]}^T.$

Furthermore, suppose $\widehat{m}\left(x\right)$ is the computed cubic B-spline approximation to $\tilde{m}\left(x\right)$ given by

(22)$\widehat{w}\left(x_i\right)=\widehat{m}\left(x_i\right)=\sum _{i=-1}^{N+1}\widehat{{\sigma }_i}{B}_i{\left(x\right)}_i$

where

(23)$\widehat{\sigma }=\widehat{{\sigma }_i}={[{\widehat{\sigma }}_{-1},{\widehat{\sigma }}_0,\dots ,{\widehat{\sigma }}_{N+1}]}^T.$

To approximate the error

$\parallel \tilde{w}\left(x_i\right)-\tilde{m}\left(x_i\right){\parallel }_{\infty },$

we need to estimate the error

${\parallel \tilde{w}\left(x_i\right)-\widehat{m}\left(x_i\right)\parallel }_{\infty }$

and

${\parallel \widehat{w}\left(x_i\right)-\tilde{m}\left(x_i\right)\parallel }_{\infty }$

differently. We know that the system of the $(n+3)\times (n+3)$ matrix leads to

(24)$A\sigma =F.$

It follows that

(25)$A\tilde{\sigma }=\tilde{F}$

and

(26)$A\widehat{\sigma }=\widehat{F}.$

Then, by substracting (25) and (26), we have

(27)$A(\widehat{\sigma }-\tilde{\sigma })=\widehat{F}-\tilde{F},$

where A is an $(n+3)\times (n+3)$ dimensional matrix, and

(28)$F={[F_{-1},F_0,\dots ,F_{N+1}]}^T,$

where T stands for transpose. Hence, from (27), we have

(29)$(\widehat{\sigma }-\tilde{\sigma })=A^{-1}(\widehat{F}-\tilde{F}).$

Now, consider taking the infinity norm from (29), and we have

$\parallel \widehat{\sigma }-\tilde{\sigma }\parallel ={\parallel A^{-1}\parallel }_{\infty }{\parallel \widehat{F}-\tilde{F}\parallel }_{\infty }.$

Note that the B-spline basis $B_{-1},B_0,B_1,\dots ,B_{N+1}$ satisfies the following inequality

(30)$\left|\sum _{i=-1}^{N+1}\widehat{{\sigma }_i}{B}_i{\left(x\right)}_i\right|\le 1.$

Adopted from [51,52,53], we have

(31)$\parallel A^{-1}{\parallel }_{\infty }{\parallel \widehat{F}-\tilde{F}\parallel }_{\infty }\le \varphi h^2,$

(32)$\parallel \widehat{\sigma }-\tilde{\sigma }{\parallel }_{\infty }\le \varphi h^2.$

Additionally,

(33)$\widehat{m}\left(x_i\right)-\tilde{m}\left(x_i\right)=(\widehat{\sigma }-\tilde{\sigma })\sum _{i=-1}^{N+1}{B}_i{\left(x\right)}_i,$

(34)${∥\widehat{m}\left(x_i\right)-\tilde{m}(x_i)∥}_{\infty }={∥(\widehat{\sigma }-\tilde{\sigma })\sum _{i=-1}^{N+1}{B}_i{\left(x\right)}_i∥}_{\infty }.$

Now, consider

(35)${∥\widehat{m}\left(x_i\right)-\tilde{m}(x_i)∥}_{\infty }\le {∥(\widehat{\sigma }-\tilde{\sigma })∥}_{\infty }\left|\sum _{i=-1}^{N+1}{B}_i{\left(x\right)}_i\right|\le \varphi h^2,$

(36)${∥\tilde{w}\left(x_i\right)-\widehat{m}(x_i)∥}_{\infty }\le \rho h^4$

and

(37)${∥\tilde{w}\left(x_i\right)-\tilde{m}(x_i)∥}_{\infty }\le {∥\tilde{w}\left(x_i\right)-\widehat{m}\left(x_i\right)∥}_{\infty }+{∥\widehat{m}\left(x_i\right)-\tilde{m}\left(x_i\right)∥}_{\infty }.$

Substituting (35) and (36) into (37), we have

(38)${∥\tilde{w}\left(x_i\right)-\tilde{m}(x_i)∥}_{\infty }\le \varphi h^2+\rho h^4=\alpha h^2,$

where $\alpha =\varphi +\rho h^2$. □

Thus, this method is second-order convergent, given by

(39)${∥\tilde{w}\left(x_i\right)-\tilde{m}(x_i)∥}_{\infty }\le \alpha h^2.$

5. Numerical Examples

In this section, three numerical problems of the linear system of ordinary differential equations are compared with the exact solutions and existing methods to demonstrate the efficiency and accuracy of the proposed method. The numerical errors are measured using the error norm $L_{\infty }$, defined as follows:

$L_{\infty }=\underset{j}{max}\left|W\left(x_j\right)-w\left(x_j\right)\right|\mathrm{or}{L}_{\infty }=\underset{j}{max}\left|Z\left(x_j\right)-z\left(x_j\right)\right|.$

5.1. Problem 1

Consider the system of the linear two-point BVPs equation [40]

(40)$\begin{array}{cc}\hfill & w^{″}\left(x\right)+xw\left(x\right)+xz\left(x\right)=2,\hfill \\ \hfill & z^{″}\left(x\right)+2xz\left(x\right)+2xw\left(x\right)=-2,\hfill \end{array}$

with boundary conditions

$\begin{array}{c}\hfill w\left(0\right)=w\left(1\right)=z\left(0\right)=0=z\left(1\right)=0.\end{array}$

The exact solutions are $w\left(x\right)=x^2-x$ and $z\left(x\right)=x-x^2$.

In Table 2, the approximate solutions, the exact solutions and the absolute errors for Problem 1 when $N=5$ are reported. It clearly shows that the approximate solutions are in good agreement with the exact solution. Table 3 lists the comparison of error norm between ECBM and the proposed method for Problem 1 with different step sizes, $N=5$ and $N=21$. Two free parameters that are involved in ECBM for this case in Table 3 are obtained by trial and error. The truncation error of the proposed method was $O\left(h^5\right)$ accurate. On the other hand, CBM and ECBM in [40] were $O\left(h^2\right)$ accurate. Our presented method produced more accurate results compared to the earlier methods.

5.2. Problem 2

Consider the system of the linear two-point BVPs equation [40]

(41)$\begin{array}{cc}\hfill & w^{″}\left(x\right)+w^{\prime }\left(x\right)+xw\left(x\right)+z^{\prime }\left(x\right)+2xz\left(x\right)=f_1\left(x\right),\hfill \\ \hfill & z^{″}\left(x\right)+z\left(x\right)+2w^{\prime }\left(x\right)+x^{2}w\left(x\right)=f_2\left(x\right),\hfill \end{array}$

with boundary conditions

$\begin{array}{c}\hfill w\left(0\right)=w\left(1\right)=z\left(0\right)=0=z\left(1\right)=0,\end{array}$

where $x\in [0,1]$, $f_1\left(x\right)=-2(x+1)cos\left(x\right)+\pi cos\left(\pi x\right)+2xsin\left(\pi x\right)+(4x-2x^2-4)sin\left(x\right)$ and $f_2\left(x\right)=-4(x-1)cos\left(x\right)-2(2-x^2+x^3)sin\left(x\right)-({\pi }^2-1)sin\left(\pi x\right)$. The true solutions are $w\left(x\right)=2(1-x)sin\left(x\right)$ and $z\left(x\right)=sin\left(\pi x\right)$.

Table 4 lists the approximate solution, the exact solution and the absolute errors for Problem 2 when $N=5$. It clearly shows that the approximate solutions promise a good agreement with the exact solution.

The comparison of the absolute errors between the present methods and the proposed method are shown for Problem 2 in Table 5 and Table 6. The reproducing kernel method in [19] was solved using eleven points in $[0,1]$, while the sinc-collocation method in [17] was solved for $N=5$ with the same number of points in $[0,1]$. The error bounds for the reproducing kernel method are at least $O\left(h^2\right)$ and $O\left(h\right)$. Two trial and error-free parameters involved in ECBM for this case are ${\lambda }_1=-1.0\times {10}^{-3}$ and ${\lambda }_2=-1.0\times {10}^{-3}$, respectively. Evidently, our proposed method produced better approximations compared with the earlier methods.

In Table 7 and Table 8, the comparison of error norm between ECBM and the presented method are tabulated for Problem 2 when $N=5$ and $N=25$, respectively. This clearly shows our presented method is more powerful. Two free parameters that are involved in ECBM for this case in Table 7 and Table 8 are obtained from the optimization technique. Table 9 reports the $L_{\infty }$ error norm with different N for Problem 2.

5.3. Problem 3

Consider the system of the linear two-point BVPs equation [21]

(42)$\begin{array}{cc}\hfill & w^{″}\left(x\right)+(2x-1)w^{\prime }\left(x\right)+cos\left(\pi x\right)z^{\prime }\left(x\right)=f_1\left(x\right),\hfill \\ \hfill & z^{″}\left(x\right)+xw\left(x\right)=f_2\left(x\right),\hfill \end{array}$

with boundary conditions

$\begin{array}{c}\hfill w\left(0\right)=w\left(1\right)=z\left(0\right)=z\left(1\right)=0,\end{array}$

where $x\in (0,1)$, $f_1\left(x\right)=-{\pi }^{2}sin\left(\pi x\right)+(2x-1)\pi cos\left(\pi x\right)+(2x-1)cos\left(\pi x\right)$ and $f_2\left(x\right)=2+xsin\left(\pi x\right).$ The analytical solutions are $w\left(x\right)=sin\left(\pi x\right)$ and $z\left(x\right)=x^2-x.$

Table 10 lists the approximate, the exact solution and the absolute error when $N=5$. It clearly shows that the approximate solutions exhibit a good agreement with the exact solution.

The comparison of the absolute errors between the existing methods and the proposed method is shown for Problem 3 when $N=41$ and $N=20$ in Table 11 and Table 12. From both tables, we noted that for $N=20$, NCBM can already match the accuracy of the VIM, CBM and ECBM. Two free parameters involved in ECBM in Table 11 and Table 12 are obtained from trial and error. The VIM in [21] was solved using one iteration step. It is observed that our proposed method is more precise compared to all earlier methods.

Table 13 compares the error norm between He’s Homotopy Perturbation, Laplace Homotopy, ECBM and the presented method for Problem 3 when $N=5$, while Table 14 compares the error norm between ECBM and the presented method for $N=41$. Two free parameters that are involved in ECBM in Table 13 are obtained from the optimization technique. Conversely, trial and error are applied to find two free parameters in Table 14. The $L_{\infty }$ error norm with different N for Problem 3 are tabulated in Table 15. It is noted that our presented method is found to be reasonably good.

Table 2, Table 4, Table 5 and Table 6 and Table 10, Table 11 and Table 12 list the numerical results, the exact solution and the absolute errors for each problem at uniform mesh. Table 9 and Table 15 present the $L_{\infty }$ norm with different N for Problem 2 and Problem 3, respectively. Additionally, for each problem, the details of error norm $L_{\infty }$ of the existing and the proposed methods at different values of N are reported in Table 3, Table 7, Table 8, Table 13 and Table 14. Consequently, the approximation obtained by the proposed method was more precise compared to others. A larger N offers greater precision but at the cost of a longer computation. This method also does not require a free parameter, but it is still the most superior and reliable method compared to the stated existing methods.

6. Conclusions

The NCBM has been applied and analyzed to numerically solve a linear system of two-point boundary value problems in this study. The method presented was based on a typical cubic B-spline, a CBS basis function that engages with the new approximation for the second-order derivative. Theoretically, it has been discovered that our method is second-order convergence. Three numerical examples were presented, and error norms, $L_{\infty }$, were calculated. We found that as the step size decreased, the error decreased, resulting in higher accuracy. Thus, it is concluded that our method gives comparable results to the stated existing methods. This method has the following three advantages: (a) it can avoid the unnecessary calculation in finding the unknown parameter; (b) it can produce up to a fifth-order accuracy; and (c) it can solve the linear system of a two-point boundary value problem accurately and efficiently. In the future, the proposed method can be used to solve more difficult problems in engineering and sciences through a graphics processing unit, GPU.

Author Contributions

Conceptualization, B.L., S.A.A.K. and I.H.; methodology, B.L., M.Y.M., S.A.A.K. and I.H.; software, B.L., M.Y.M. and S.A.A.K.; validation, B.L., M.Y.M., S.A.A.K. and I.H.; formal analysis, B.L., M.Y.M., S.A.A.K. and I.H.; investigation, B.L.; writing the original draft preparation, B.L.; writing the review and editing, B.L., S.A.A.K. and I.H.; supervision, S.A.A.K. and I.H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The authors confirm that the data supporting the findings of this study are available within the article.

Acknowledgments

The authors gratefully acknowledge the anonymous referees for their insightful comments and suggestions that helped to enhance this paper. B. Latif is thankful for the full scholarship (SLAB) awarded by the Malaysian Ministry of Higher Education (MOHE) and Universiti Teknologi MARA (UiTM) Malaysia.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

BVPs	Boundary value-problems
NCBM	New symmetric cubic B-spline method
IVPs	Initial value-problems
CBS	Cubic B-spline
CBM	Cubic B-spline method
ECBM	Extended cubic B-spline method
VIM	Variational iteration method

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Tables

Table 1

Coefficient of $B_i\left(x\right),B_i^{\prime }\left(x\right)$ and $B_i^{″}\left(x\right)$ at the nodes.

	${\mathit{x}}_{\mathit{i}-2}$	${\mathit{x}}_{\mathit{i}-1}$	${\mathit{x}}_{\mathit{i}}$	${\mathit{x}}_{\mathit{i}+1}$	${\mathit{x}}_{\mathit{i}+2}$
$B_i\left(x\right)$	0	1	4	1	0
$B_i^{\prime }\left(x\right)$	0	$-1/2h$	0	$1/2h$	0
$B_i^{″}\left(x\right)$	0	$1/h^2$	$-2/h^2$	$1/h^2$	0

Table 2

Absolute errors for Problem 1 when $N=5$.

x	NCBM	Exact	Absolute Error	NCBM	Exact	Absolute Error
	$\mathit{W}\left(\mathit{x}\right)$	Solution $\mathit{w}\left(\mathit{x}\right)$	$\mathit{W}\left(\mathit{x}\right)-\mathit{w}\left(\mathit{x}\right)$	$\mathit{Z}\left(\mathit{x}\right)$	Solution $\mathit{z}\left(\mathit{x}\right)$	$\mathit{Z}\left(\mathit{x}\right)-\mathit{z}\left(\mathit{x}\right)$
0.2	$-0.16$	$-0.16$	$2.78\times {10}^{-17}$	0.16	0.16	$2.22\times {10}^{-16}$
0.4	$-0.24$	$-0.24$	$1.11\times {10}^{-16}$	0.24	0.24	$5.55\times {10}^{-17}$
0.6	$-0.24$	$-0.24$	$1.67\times {10}^{-16}$	0.24	0.24	$1.11\times {10}^{-16}$
0.8	$-0.16$	$-0.16$	$1.39\times {10}^{-16}$	0.16	0.16	$2.78\times {10}^{-17}$

Table 3

The $L_{\infty }$ error norm for Problem 1 when $N=5$ and $N=21$.

	ECBM [40]	NCBM	ECBM [40]	ECBM [40]	NCBM
	$\mathit{N}=\mathbf{5}$	$\mathit{N}=\mathbf{5}$	$\mathit{N}=\mathbf{21}$	$\mathit{N}=\mathbf{21}$	$\mathit{N}=\mathbf{21}$
	${\mathbf{\lambda }}_{\mathbf{1}}={\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{0}$		${\mathbf{\lambda }}_{\mathbf{1}}={\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{0}$	${\mathbf{\lambda }}_{\mathbf{1}}={\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{1.25}\times {\mathbf{10}}^{-\mathbf{14}}$
$W\left(x\right)$	$3.47\times {10}^{-15}$	$1.67\times {10}^{-16}$	$3.72\times {10}^{-13}$	$1.73\times {10}^{-13}$	$1.07\times {10}^{-15}$
$Z\left(x\right)$	$3.69\times {10}^{-15}$	$2.22\times {10}^{-16}$	$2.53\times {10}^{-13}$	$1.67\times {10}^{-13}$	$6.94\times {10}^{-16}$

Table 4

Absolute errors for Problem 2 when $N=5$.

x	NCBM	Exact	Absolute Error	NCBM	Exact	Absolute Error
	$\mathit{W}\left(\mathit{x}\right)$	Solution $\mathit{w}\left(\mathit{x}\right)$	$\mathit{W}\left(\mathit{x}\right)-\mathit{w}\left(\mathit{x}\right)$	$\mathit{Z}\left(\mathit{x}\right)$	Solution $\mathit{z}\left(\mathit{x}\right)$	$\mathit{Z}\left(\mathit{x}\right)-\mathit{z}\left(\mathit{x}\right)$
0.2	0.317794	0.317871	$7.71\times {10}^{-5}$	0.587778	0.587785	$6.92\times {10}^{-6}$
0.4	0.467274	0.467302	$2.81\times {10}^{-5}$	0.950901	0.951057	$1.55\times {10}^{-4}$
0.6	0.451776	0.451714	$6.19\times {10}^{-5}$	0.950901	0.951057	$1.56\times {10}^{-4}$
0.8	0.287036	0.286942	$9.37\times {10}^{-5}$	0.587778	0.587785	$6.98\times {10}^{-6}$

Table 5

Absolute errors for Problem 2 for $w\left(x\right)$.

x			ECBM [40]	ECBM [40]
	Reproducing	Sinc−	${\mathit{\lambda }}_1=0$	${\mathit{\lambda }}_1=-1.0\times {10}^{-3}$	NCBM
	Kernel [19]	Collocation [17]	${\mathit{\lambda }}_2=0$	${\mathit{\lambda }}_2=-1.0\times {10}^{-3}$	$\mathit{N}=25$
		$\mathit{N}=\mathbf{5}$	$\mathit{N}=\mathbf{25}$	$\mathit{N}=\mathbf{25}$
0.08	$3.3\times {10}^{-3}$	$3.2\times {10}^{-3}$	$1.3\times {10}^{-4}$	$1.4\times {10}^{-5}$	$9.6\times {10}^{-8}$
0.24	$7.7\times {10}^{-3}$	$9.2\times {10}^{-4}$	$2.7\times {10}^{-4}$	$1.1\times {10}^{-5}$	$1.6\times {10}^{-7}$
0.4	$9.7\times {10}^{-3}$	$2.0\times {10}^{-3}$	$2.7\times {10}^{-4}$	$2.1\times {10}^{-5}$	$9.2\times {10}^{-8}$
0.56	$9.5\times {10}^{-3}$	$2.2\times {10}^{-4}$	$2.0\times {10}^{-4}$	$5.9\times {10}^{-5}$	$1.9\times {10}^{-8}$
0.72	$7.3\times {10}^{-3}$	$4.1\times {10}^{-3}$	$9.4\times {10}^{-5}$	$7.8\times {10}^{-5}$	$1.0\times {10}^{-7}$
0.88	$3.4\times {10}^{-3}$	$1.0\times {10}^{-2}$	$1.6\times {10}^{-5}$	$5.6\times {10}^{-5}$	$9.3\times {10}^{-8}$
0.96	$1.1\times {10}^{-3}$	$2.1\times {10}^{-3}$	$3.6\times {10}^{-8}$	$2.3\times {10}^{-5}$	$4.0\times {10}^{-8}$

Table 6

Absolute errors for Problem 2 for $z\left(x\right)$.

x			ECBM [40]	ECBM [40]
	Reproducing	Sinc−	${\mathit{\lambda }}_1=0$	${\mathit{\lambda }}_1=-1.0\times {10}^{-3}$	NCBM
	Kernel [19]	Collocation [17]	${\mathit{\lambda }}_2=0$	${\mathit{\lambda }}_2=-1.0\times {10}^{-3}$	$\mathit{N}=25$
		$\mathit{N}=\mathbf{5}$	$\mathit{N}=\mathbf{25}$	$\mathit{N}=\mathbf{25}$
0.08	$7.7\times {10}^{-3}$	$1.5\times {10}^{-3}$	$3.8\times {10}^{-4}$	$2.2\times {10}^{-4}$	$1.8\times {10}^{-7}$
0.24	$2.0\times {10}^{-2}$	$7.0\times {10}^{-3}$	$9.9\times {10}^{-4}$	$6.0\times {10}^{-4}$	$4.7\times {10}^{-7}$
0.4	$2.7\times {10}^{-2}$	$7.4\times {10}^{-3}$	$1.3\times {10}^{-3}$	$8.3\times {10}^{-4}$	$6.3\times {10}^{-7}$
0.56	$2.7\times {10}^{-2}$	$1.0\times {10}^{-2}$	$1.4\times {10}^{-3}$	$8.6\times {10}^{-4}$	$6.4\times {10}^{-7}$
0.72	$2.0\times {10}^{-2}$	$4.4\times {10}^{-3}$	$1.1\times {10}^{-3}$	$6.8\times {10}^{-4}$	$5.1\times {10}^{-7}$
0.88	$9.4\times {10}^{-3}$	$2.1\times {10}^{-2}$	$5.0\times {10}^{-4}$	$3.3\times {10}^{-4}$	$2.5\times {10}^{-7}$
0.96	$3.1\times {10}^{-3}$	$6.9\times {10}^{-3}$	$1.7\times {10}^{-4}$	$1.1\times {10}^{-4}$	$8.5\times {10}^{-8}$

Table 7

The $L_{\infty }$ error norm for Problem 2 when $N=5$.

	ECBM [40]	ECBM [40]	NCBM
	${\mathbf{\lambda }}_{\mathbf{1}}={\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{0}$	${\mathbf{\lambda }}_{\mathbf{1}}=-\mathbf{1.269208}\times {\mathbf{10}}^{-\mathbf{2}}$
		${\mathbf{\lambda }}_{\mathbf{2}}=-\mathbf{6.634523}\times {\mathbf{10}}^{-\mathbf{2}}$
$W\left(x\right)$	$2.09\times {10}^{-3}$	$1.80\times {10}^{-5}$	$9.37\times {10}^{-5}$
$Z\left(x\right)$	$1.75\times {10}^{-4}$	$1.75\times {10}^{-4}$	$1.56\times {10}^{-4}$

Table 8

The $L_{\infty }$ error norm for Problem 2 when $N=25$.

	ECBM [40]	ECBM [40]	NCBM
	${\mathbf{\lambda }}_{\mathbf{1}}={\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{0}$	${\mathbf{\lambda }}_{\mathbf{1}}=-\mathbf{1.0}\times {\mathbf{10}}^{-\mathbf{3}}$
		${\mathbf{\lambda }}_{\mathbf{2}}=-\mathbf{1.0}\times {\mathbf{10}}^{-\mathbf{3}}$
$W\left(x\right)$	$2.72\times {10}^{-4}$	$7.80\times {10}^{-5}$	$1.56\times {10}^{-7}$
$Z\left(x\right)$	$1.36\times {10}^{-3}$	$8.60\times {10}^{-4}$	$6.53\times {10}^{-7}$

Table 9

The $L_{\infty }$ error norm with different N for Problem 2.

N	NCBM $\mathit{W}\mathbf{\left(}\mathit{x}\mathbf{\right)}$	NCBM $\mathit{Z}\mathbf{\left(}\mathit{x}\mathbf{\right)}$
40	$2.40\times {10}^{-8}$	$1.00\times {10}^{-7}$
80	$1.50\times {10}^{-9}$	$6.27\times {10}^{-9}$
100	$6.14\times {10}^{-10}$	$2.57\times {10}^{-9}$

Table 10

Absolute errors for Problem 3 when $N=5$.

x	NCBM	Exact	Absolute Error	NCBM	Exact	Absolute Error
	$\mathit{W}\left(\mathit{x}\right)$	Solution $\mathit{w}\left(\mathit{x}\right)$	$\mathit{W}\left(\mathit{x}\right)-\mathit{w}\left(\mathit{x}\right)$	$\mathit{Z}\left(\mathit{x}\right)$	Solution $\mathit{z}\left(\mathit{x}\right)$	$\mathit{Z}\left(\mathit{x}\right)-\mathit{z}\left(\mathit{x}\right)$
0.2	0.587845	0.587785	$5.98\times {10}^{-5}$	$-0.160001$	$-0.160000$	$1.10\times {10}^{-6}$
0.4	0.950959	0.951057	$9.72\times {10}^{-5}$	$-0.240002$	$-0.240000$	$2.47\times {10}^{-6}$
0.6	0.950959	0.951057	$9.72\times {10}^{-5}$	$-0.240002$	$-0.240000$	$2.37\times {10}^{-6}$
0.8	0.587845	0.587785	$5.97\times {10}^{-5}$	$-0.160000$	$-0.160000$	$4.86\times {10}^{-7}$

Table 11

Absolute errors for Problem 3 for $w\left(x\right)$.

x				ECBM [40]
	VIM	CBM	ECBM [40]	${\mathit{\lambda }}_1=-1.0\times {10}^{-3}$	NCBM	NCBM
	[21]	[40]	${\mathit{\lambda }}_1=0$	${\mathit{\lambda }}_2=0$	$\mathit{N}=41$	$\mathit{N}=20$
		$\mathit{N}=\mathbf{41}$	${\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{0}$	$\mathit{N}=\mathbf{41}$
0.1	$3.30\times {10}^{-4}$	$1.40\times {10}^{-4}$	$1.30\times {10}^{-4}$	$2.83\times {10}^{-6}$	$1.63\times {10}^{-8}$	$2.83\times {10}^{-7}$
0.2	$2.51\times {10}^{-3}$	$2.80\times {10}^{-4}$	$2.56\times {10}^{-4}$	$5.55\times {10}^{-6}$	$3.57\times {10}^{-8}$	$6.34\times {10}^{-7}$
0.3	$7.84\times {10}^{-3}$	$3.90\times {10}^{-4}$	$3.60\times {10}^{-4}$	$7.81\times {10}^{-6}$	$5.39\times {10}^{-8}$	$9.59\times {10}^{-7}$
0.4	$1.66\times {10}^{-2}$	$4.60\times {10}^{-4}$	$4.28\times {10}^{-4}$	$9.30\times {10}^{-6}$	$6.71\times {10}^{-8}$	$1.19\times {10}^{-6}$
0.5	$2.77\times {10}^{-2}$	$4.80\times {10}^{-2}$	$4.52\times {10}^{-4}$	$9.82\times {10}^{-6}$	$7.26\times {10}^{-8}$	$1.27\times {10}^{-6}$
0.6	$3.87\times {10}^{-2}$	$4.60\times {10}^{-2}$	$4.28\times {10}^{-4}$	$9.30\times {10}^{-6}$	$6.93\times {10}^{-8}$	$1.19\times {10}^{-6}$
0.7	$4.59\times {10}^{-2}$	$3.90\times {10}^{-2}$	$3.60\times {10}^{-4}$	$7.81\times {10}^{-6}$	$5.78\times {10}^{-8}$	$9.59\times {10}^{-7}$
0.8	$4.49\times {10}^{-2}$	$2.80\times {10}^{-2}$	$2.56\times {10}^{-4}$	$5.56\times {10}^{-6}$	$4.05\times {10}^{-8}$	$6.34\times {10}^{-7}$
0.9	$3.09\times {10}^{-2}$	$1.50\times {10}^{-2}$	$1.30\times {10}^{-4}$	$2.83\times {10}^{-6}$	$2.10\times {10}^{-8}$	$6.83\times {10}^{-7}$

Table 12

Absolute errors for Problem 3 for $z\left(x\right)$.

x		ECBM [40]	ECBM [40]
	CBM	${\mathit{\lambda }}_1=0$	${\mathit{\lambda }}_1=-1.0\times {10}^{-3}$	NCBM	NCBM
	[40]	${\mathit{\lambda }}_2=0$	${\mathit{\lambda }}_2=0$	$\mathit{N}=41$	$\mathit{N}=20$
	$\mathit{N}=\mathbf{41}$	$\mathit{N}=\mathbf{41}$	$\mathit{N}=\mathbf{41}$
0.1	$5.74\times {10}^{-6}$	$5.74\times {10}^{-6}$	$1.25\times {10}^{-7}$	$8.64\times {10}^{-10}$	$1.54\times {10}^{-8}$
0.2	$1.13\times {10}^{-5}$	$1.13\times {10}^{-5}$	$2.46\times {10}^{-7}$	$1.71\times {10}^{-9}$	$3.05\times {10}^{-8}$
0.3	$1.64\times {10}^{-5}$	$1.64\times {10}^{-5}$	$3.56\times {10}^{-7}$	$2.49\times {10}^{-9}$	$4.42\times {10}^{-8}$
0.4	$2.03\times {10}^{-5}$	$2.03\times {10}^{-5}$	$4.42\times {10}^{-7}$	$3.11\times {10}^{-9}$	$5.51\times {10}^{-8}$
0.5	$2.26\times {10}^{-5}$	$2.26\times {10}^{-5}$	$4.91\times {10}^{-7}$	$3.49\times {10}^{-9}$	$6.12\times {10}^{-8}$
0.6	$2.26\times {10}^{-5}$	$2.26\times {10}^{-5}$	$4.92\times {10}^{-7}$	$3.53\times {10}^{-9}$	$6.10\times {10}^{-8}$
0.7	$2.01\times {10}^{-5}$	$2.01\times {10}^{-5}$	$4.37\times {10}^{-7}$	$3.20\times {10}^{-9}$	$5.39\times {10}^{-8}$
0.8	$1.51\times {10}^{-5}$	$1.51\times {10}^{-5}$	$3.29\times {10}^{-7}$	$2.49\times {10}^{-9}$	$4.01\times {10}^{-8}$
0.9	$8.14\times {10}^{-6}$	$8.14\times {10}^{-6}$	$1.76\times {10}^{-7}$	$1.48\times {10}^{-9}$	$2.13\times {10}^{-8}$

Table 13

The $L_{\infty }$ error norm for Problem 3 when $N=5$.

			ECBM [40]	ECBM [40]
	He’s Homotopy	Laplace	${\mathit{\lambda }}_1=0$	${\mathit{\lambda }}_1=-6.639145\times {10}^{-2}$	NCBM
	Pertubation [18]	Homotopy [13]	${\mathit{\lambda }}_2=0$	${\mathit{\lambda }}_2=1.161882\times {10}^{-6}$
$W\left(x\right)$	$2.1\times {10}^{-4}$	$2.2\times {10}^{-5}$	$2.8\times {10}^{-2}$	$1.4\times {10}^{-4}$	$9.7\times {10}^{-5}$
$Z\left(x\right)$	$3.2\times {10}^{-4}$	$1.1\times {10}^{-5}$	$1.4\times {10}^{-3}$	$7.2\times {10}^{-6}$	$2.5\times {10}^{-6}$

Table 14

The $L_{\infty }$ error norm for Problem 3 when $N=41$.

	ECBM [40]	ECBM [40]
	${\mathbf{\lambda }}_{\mathbf{1}}=\mathbf{0}$	${\mathbf{\lambda }}_{\mathbf{1}}=-\mathbf{1.0}\times {\mathbf{10}}^{-\mathbf{3}}$	NCBM
	${\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{0}$	${\mathbf{\lambda }}_{\mathbf{2}}=\mathbf{0}$
$W\left(x\right)$	$4.52\times {10}^{-4}$	$9.82\times {10}^{-6}$	$7.26\times {10}^{-8}$
$Z\left(x\right)$	$2.26\times {10}^{-5}$	$4.92\times {10}^{-7}$	$3.56\times {10}^{-9}$

Table 15

The $L_{\infty }$ error norm with different N for Problem 3.

N	NCBM	NCBM
	$\mathit{W}\mathbf{\left(}\mathit{x}\mathbf{\right)}$	$\mathit{Z}\mathbf{\left(}\mathit{x}\mathbf{\right)}$
60	$1.59\times {10}^{-8}$	$7.77\times {10}^{-10}$
80	$5.02\times {10}^{-9}$	$2.46\times {10}^{-10}$
100	$2.06\times {10}^{-9}$	$1.01\times {10}^{-10}$

References

1. Ha, S.N. A nonlinear shooting method for two-point boundary value problems. Comput. Math. Appl.; 2001; 42, pp. 1411-1420. [DOI: https://dx.doi.org/10.1016/S0898-1221(01)00250-4]

2. Chawla, M.M. An eighth order tridiagonal finite difference method for nonlinear two-point boundary value problems. Bit Numer. Math.; 1977; 17, pp. 281-285. [DOI: https://dx.doi.org/10.1007/BF01932148]

3. Usmani, R.A. A method of high-order accuracy for the numerical integration of boundary value problems. BIT; 1973; 13, pp. 458-469. [DOI: https://dx.doi.org/10.1007/BF01933410]

4. Tirmizi, I.A.; Twizelle, E.H. Higher-order finite-difference methods for nonlinear second-order two-point boundary-value problems. Appl. Math. Lett.; 1977; 15, pp. 897-902. [DOI: https://dx.doi.org/10.1016/S0893-9659(02)00060-5]

5. Boyd, J.P. Chebyshev and Fourier Spectral Methods; 2nd ed. Courier Corporation: Mineola, NY, USA, 2000.

6. Mehrpouya, M.A.; Shamsi, M. Gauss pseudospectral and continuation methods for solving two-point boundary value problems in optimal control theory. Appl. Math. Model.; 2015; 39, pp. 5047-5057. [DOI: https://dx.doi.org/10.1016/j.apm.2015.04.009]

7. Gheorghiu, C.I. Chebfun solutions to a class of 1D singular and nonlinear boundary value problems. Computation; 2022; 10, 116. [DOI: https://dx.doi.org/10.3390/computation10070116]

8. Betts, J.T. Survey of numerical methods for trajectory optimization. J. Guid. Control. Dyn.; 1998; 21, pp. 193-207. [DOI: https://dx.doi.org/10.2514/2.4231]

9. Wang, H.Z.; Li, Y. Existence and uniqueness of solutions to two-point boundary value problems for ordinary differential equations. Z. Angew. Math. Phys. ZAMP; 1996; 47, pp. 373-384. [DOI: https://dx.doi.org/10.1007/BF00916644]

10. Nikooeinejad, Z.; Heydari, M.; Loghmani, G.B. A numerical iterative method for solving two-point BVPs in infinite-horizon nonzero–sum differential games: Economic applications. Math. Comput. Simul.; 2022; 200, pp. 404-427. [DOI: https://dx.doi.org/10.1016/j.matcom.2022.04.022]

11. Sweidan, M.; Chen, X.; Zheng, X. The Shortley–Weller scheme for variable coefficient two-point boundary value problems and its application to tumor growth problem with heterogeneous microenvironment. J. Comput. Appl. Math.; 2020; 376, 112874. [DOI: https://dx.doi.org/10.1016/j.cam.2020.112874]

12. Campos, L.M.B.C.; Gil, P.J.S. The two-point boundary-value problem for rocket trajectories. Aerospace; 2020; 7, 131. [DOI: https://dx.doi.org/10.3390/aerospace7090131]

13. Ogunlaran, O.M.; Ademola, A.T. On the Laplace homotopy analysis method for a non-linear system of second-order boundary value problems. Gen. Math. Notes; 2015; 26, pp. 11-22.

14. Sami Bataineh, A.; Noorani, M.S.M.; Hashim, I. Modified homotopy analysis method for solving systems of second-order BVPs. Commun. Nonlinear Sci. Numer. Simul.; 2009; 14, pp. 430-442. [DOI: https://dx.doi.org/10.1016/j.cnsns.2007.09.012]

15. Arqub, O.A.; Abo-Hammour, Z. Numerical solution of systems of second-order boundary value problems using continuous genetic algorithm. Inf. Sci.; 2014; 279, pp. 396-415. [DOI: https://dx.doi.org/10.1016/j.ins.2014.03.128]

16. El-Gamel, M. Sinc-collocation method for solving linear and non-linear system of second-order boundary value problems. Appl. Math.; 2012; 3, pp. 1627-1633. [DOI: https://dx.doi.org/10.4236/am.2012.311225]

17. Dehghan, M.; Saadatmandi, A. The numerical solution of a non-linear system of second-order boundary value problems using the sinc-collocation method. Math. Comput. Model.; 2007; 46, pp. 1434-1441. [DOI: https://dx.doi.org/10.1016/j.mcm.2007.02.002]

18. Saadatmandi, A.; Dehghan, M.; Eftekhari, A. Application of H e’s homotopy perturbation method for non-linear system of second-order boundary value problems. Nonlinear Anal. Real World Appl.; 2009; 10, pp. 1912-1922. [DOI: https://dx.doi.org/10.1016/j.nonrwa.2008.02.032]

19. Geng, F.; Cui, M. Solving a non-linear system of second order boundary value problems. J. Math. Anal. Appl.; 2007; 327, pp. 1167-1181. [DOI: https://dx.doi.org/10.1016/j.jmaa.2006.05.011]

20. Ratib Anakira, N.; Alomari, A.K.; Jameel, A.F.; Hashim, I. Multistage optimal homotopy asymptotic method for solving initial–value problems. J. Nonlinear Sci. Appl.; 2016; 9, pp. 1826-1843. [DOI: https://dx.doi.org/10.22436/jnsa.009.04.37]

21. Lu, J. Variational iteration method for solving a non-linear system of second-order boundary value problems. Comput. Math. Appl.; 2007; 54, pp. 1133-1138. [DOI: https://dx.doi.org/10.1016/j.camwa.2006.12.060]

22. Saadatmandi, A.; Askari Farsangi, J. Chebyshev finite difference method for a non-linear system of second-order boundary value problems. Appl. Math. Comput.; 2007; 192, pp. 586-591.

23. Akram, T.; Abbas, M.; Iqbal, A.; Baleanu, D.; Asad, H.H. Novel numerical approach based on modified extended cubic B–spline functions for solving non-linear time-fractional telegraph equation. Symmetry; 2020; 12, 1154. [DOI: https://dx.doi.org/10.3390/sym12071154]

24. Tayebi, S.; Momani, S.; Abu Arqub, O. The cubic B-spline interpolation method for numerical point solutions of conformable boundary value problems. Alex. Eng. J.; 2022; 61, pp. 1519-1528. [DOI: https://dx.doi.org/10.1016/j.aej.2021.06.057]

25. Parumasur, N.; Adetona, R.A.; Singh, P. Efficient solution of burgers’, modified burgers’ and KdV–burgers’ equations using B-spline approximation functions. Mathematics; 2023; 11, 1847. [DOI: https://dx.doi.org/10.3390/math11081847]

26. Abbas, M.; Bibi, A.; Alzaidi, A.S.M.; Nazir, T.; Majeed, A.; Akram, G. Numerical solutions of third-order time-fractional differential equations using cubic B-Spline functions. Fractal Fract.; 2022; 6, 528. [DOI: https://dx.doi.org/10.3390/fractalfract6090528]

27. Arqub, O.A.; Tayebi, S.; Baleanu, D.; Osman, M.S.; Mahmoud, W.; Alsulam, H. A numerical combined algorithm in cubic B-spline method and finite difference technique for the time-fractional nonlinear diffusion wave equation with reaction and damping terms. Results Phys.; 2022; 6, 105912. [DOI: https://dx.doi.org/10.1016/j.rinp.2022.105912]

28. Sakai, M. Piecewise cubic interpolation and two-point boundary value problems. Publ. Res. Inst. Math.; 1971; 7, pp. 345-362. [DOI: https://dx.doi.org/10.2977/prims/1195193546]

29. Albasiny, E.L.; Hoskins, W.D. Cubic spline solutions to two-point boundary value problems. Comput. J.; 1969; 12, pp. 151-153. [DOI: https://dx.doi.org/10.1093/comjnl/12.2.151]

30. Al-Said, E.A. Cubic spline method for solving two-point boundary-value problems. Korean J. Comput. Appl. Math.; 1998; 5, pp. 669-680. [DOI: https://dx.doi.org/10.1007/BF03008890]

31. Müllenheim, G. Solving two-point boundary value problems with spline functions. Ima J. Numer. Anal.; 1992; 12, pp. 503-518. [DOI: https://dx.doi.org/10.1093/imanum/12.4.503]

32. Khalifa, A.K.A.; Eilbeck, J.C. Collocation with quadratic and cubic splines. Ima J. Numer. Anal.; 1982; 2, pp. 111-121. [DOI: https://dx.doi.org/10.1093/imanum/2.1.111]

33. Al-Said, E.A.; Noor, M.A.; Al-Shejari, A.A. Numerical solutions for system of second-order boundary value problems. Korean J. Comput. Appl. Math.; 1998; 5, pp. 659-667. [DOI: https://dx.doi.org/10.1007/BF03008889]

34. Caglar, H.; Caglar, N.; Elfaituri, K. B-spline interpolation compared with finite difference, finite element and finite volume methods which applied to two-point boundary value problems. Appl. Math. Comput.; 2006; 175, pp. 72-79. [DOI: https://dx.doi.org/10.1016/j.amc.2005.07.019]

35. Hamid, N.N.A.; Majid, A.A.; Ismail, A.I.M. Extended cubic B-spline method for linear two-point boundary value problems. Sains Malays.; 2011; 40, pp. 1285-1290.

36. Hamid, N.N.A.; Majid, A.A.; Izani, A. Extended cubic B-spline interpolation method applied to linear two-point boundary value problems. World Acad. Sci. Eng. Technol.; 2010; 62, pp. 566-568.

37. Salem Heilat, A.; Ismail, A.I.M. Hybrid cubic B-spline method for solving non-linear two-point boundary value problems. Int. J. Pure Appl. Math.; 2016; 110, pp. 369-381.

38. Heilat, A.; Zureigat, H.; Batiha, B. New spline method for solving linear two-point boundary value problems. Eur. J. Pure Appl. Math.; 2021; 14, pp. 1283-1294. [DOI: https://dx.doi.org/10.29020/nybg.ejpam.v14i4.4124]

39. Caglar, N.; Caglar, H. B-spline method for solving linear system of second-order boundary value problems. Comput. Math. Appl.; 2009; 57, pp. 757-762. [DOI: https://dx.doi.org/10.1016/j.camwa.2008.09.033]

40. Heilat, A.S.; Hamid, N.A.A.; Ismail, A.I.M. Extended cubic B-spline method for solving a linear system of second-order boundary value problems. SpringerPlus; 2016; 5, pp. 1-18. [DOI: https://dx.doi.org/10.1186/s40064-016-2936-4]

41. Zhang, J.; Niu, J. Lobatto-reproducing kernel method for solving a linear system of second order boundary value problems. J. Appl. Math. Comput.; 2016; 68, pp. 3631-3653. [DOI: https://dx.doi.org/10.1007/s12190-021-01685-9]

42. Lang, F.-G.; Xu, X.-P. A new cubic B-spline method for approximating the solution of a class of non-linear second-order boundary value problem with two dependent variables. Sci. Asia; 2014; 40, pp. 444-450. [DOI: https://dx.doi.org/10.2306/scienceasia1513-1874.2014.40.444]

43. Kashif Iqbal, M.; Abbas, M.; Wasim, I. New cubic B-spline approximation for solving third order Emden–Flower type equations. Appl. Math. Comput.; 2018; 331, pp. 319-333.

44. Wasim, I.; Abbas, M.; Kashif Iqbal, M. A new extended B-spline approximation technique for second order singular boundary value problems arising in physiology. J. Math. Comput. Sci.; 2019; 19, pp. 258-267. [DOI: https://dx.doi.org/10.22436/jmcs.019.04.06]

45. Abbas, M.; Iqbal, M.K.; Zafar, B.; Zin, S.B.M. New cubic B-spline approximations for solving non-linear third-order Korteweg-De Vries equation. Indian J. Sci. Technol.; 2019; 12, pp. 1-9. [DOI: https://dx.doi.org/10.17485/ijst/2019/v12i15/141526]

46. Nazir, T.; Abbas, M.; Iqbal, M.K. A new quintic B-spline approximation for numerical treatment of Boussinesq equation. J. Math. Comput. Sci.; 2020; 20, pp. 30-42. [DOI: https://dx.doi.org/10.22436/jmcs.020.01.04]

47. Nazir, T.; Abbas, M.; Iqbal, M.K. New cubic B-spline approximation technique for numerical solutions of coupled viscous Burgers equations. Eng. Comput.; 2020; 38, pp. 83-106. [DOI: https://dx.doi.org/10.1108/EC-08-2019-0365]

48. Goh, J.; Majid, A.A.; Ismail, A.I.M. Extended cubic uniform B-spline for a class of singular boundary value problems. Scienceasia; 2011; 37, pp. 79-82. [DOI: https://dx.doi.org/10.2306/scienceasia1513-1874.2011.37.079]

49. Latif, B.; Karim, S.A.A.; Hashim, I. New cubic B-spline approximation for solving linear two-point boundary-value problems. Mathematics; 2021; 9, 1250. [DOI: https://dx.doi.org/10.3390/math9111250]

50. Tolmachev, D.; Chertovskih, R.; Zheligovsky, V. Algorithmic aspects of simulation of magnetic field generation by thermal convection in a plane layer of fluid. Mathematics; 2023; 11, 808. [DOI: https://dx.doi.org/10.3390/math11040808]

51. De Boor, C. On the convergence of odd-degree spline interpolation. J. Approx. Theory; 1968; 1, pp. 452-463. [DOI: https://dx.doi.org/10.1016/0021-9045(68)90033-6]

52. Hall, C.A. On error bounds for spline interpolation. J. Approx. Theory; 1968; 1, pp. 209-218. [DOI: https://dx.doi.org/10.1016/0021-9045(68)90025-7]

53. Prenter, P.M. Splines and Variational Methods; Courier Corporation: Dover, NY, USA, 2008.