Mathematical Problems in Engineering
Volume 2018 (2018), Article ID 2056834, 12 pages
https://doi.org/10.1155/2018/2056834

Research Article

Diagonal Block Method for Solving Two-Point Boundary Value Problems with Robin Boundary Conditions

Nadirah Mohd Nasir,1,2 Zanariah Abdul Majid,1,3 Fudziah Ismail,1,3 and Norfifah Bachok1,3

¹Institute for Mathematical Research, Universiti Putra Malaysia, 43400 Serdang, Selangor Darul Ehsan, Malaysia
²Faculty of Industrial Sciences & Technology, Universiti Malaysia Pahang, 26300 Gambang, Pahang, Malaysia
³Department of Mathematics, Faculty of Science, Universiti Putra Malaysia, 43400 Serdang, Selangor Darul Ehsan, Malaysia

Correspondence should be addressed to Zanariah Abdul Majid; zana_majid99@yahoo.com

Received 3 April 2018; Accepted 2 September 2018; Published 2 October 2018

Academic Editor: Nunzio Salerno

Copyright © 2018 Nadirah Mohd Nasir et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This numerical study presents the diagonal block method of order four for solving the second-order boundary value problems (BVPs) with Robin boundary conditions at two-point concurrently using constant step size. The solution is obtained directly without reducing to a system of first-order differential equations using a combination of predictor-corrector mode via shooting technique. The shooting method was adapted with the Newton divided difference interpolation approach as the strategy of seeking for the new initial estimate. Five numerical examples are included to examine and illustrate the practical usefulness of the proposed method. Numerical tested problem is also highlighted on the diffusion of heat generated application that imposed the Robin boundary conditions. The present findings revealed that the proposed method gives an efficient performance in terms of accuracy, total function calls, and execution time as compared with the existing method.

1. Introduction

This study is focusing on the numerical approach for solving second-order boundary value problems (BVPs) associated with Robin boundary conditions. Generally, this type of BVPs is given as follows:

with

where , and are all constants and , and all nonzero. In the case of Robin type, both functional value and derivative of the solutions are given in (2). If only functional value is given, then condition in (2) is known as Dirichlet type; otherwise BVPs will be subject to Neumann condition when only derivative values exist. Robin boundary conditions arise in several branches of applications such as in electromagnetic problem and heat transfer problem where these Robin type conditions are called impedance boundary conditions and the convective boundary conditions are, respectively, as explained in [1]. The Bernoulli polynomial together with Galerkin approximation in solving linear and nonlinear Robin boundary condition problems had been studied by Islam and Shirin [2]. The eminent scholars including Duan et al. [3] and Rach et al. [4] had derived the Adomian decomposition methods for solving BVPs imposing this condition where the approach involved analytic stage and numeric simulations. Meanwhile, Bhatta and Sastri [5] and Lang and Xu [6] had transformed respective BVPs into discretization form and solved them using symmetric global (continuous) spline and Quintic B-spline collocation method, respectively.

To optimize the computational cost, the development of proposed algorithm must at least satisfy the facility to generate solutions at several points simultaneously as suggested in Fatunla [7]. This implementation has been shown in the literature discussed by Phang et al. [8], Majid et al. [9], and Omar and Adeyeye [10]. They obtained the approximate solution of (1) at two points concurrently. Therefore, the advantages of outcome from their discussion have motivated us in this study.

Diagonal block method for solving differential equations was widely studied in the previous literatures. These include the discussion on solving first-order differential equations in [11, 12]. Solving second-order ordinary differential equations using diagonal block method has been discussed in [13]. The formulation in [12, 13] is based on backward differentiation approach. In [11], the author has derived the diagonal block method using Lagrange interpolation polynomial for solving first-order ordinary differential equations. In this research, we have extended the derivation in [11] to obtain the formulation of direct integration for solving second-order differential equations. Then the new formulae obtained have been implemented to solve the boundary value problems. The investigation will be focusing on the Newton divided difference interpolation as an iterative technique in seeking as well as updating the missing initial guesses while performing the shooting technique.

The organization of this paper is as follows. Section 2 presents the derivation of the two-point diagonal block method. Section 3 elaborates on the analysis of the method including the order, consistency, and stability. Section 4 presents the idea of the procedure of shooting method together with Newton divided difference interpolation as an iterative part. For validation and a clear overview, five tested problems will be discussed in Section 5. Finally, Section 6 concludes the finding from this study.

2. Derivation of the Diagonal Block Method

The interval of is divided into a series of blocks so that each block contains two subintervals with equal distance step size, , as depicted in Figure 1. Both numerical solutions of and in these subintervals are simultaneously generated using the appropriate number of back values. The approximate solution of at the point will be computed using three back values at the points, , and While the approximation of at the point required an additional one back value, The derivation formulae of and will be obtained by integrate equation (1) once and twice over the intervals and , respectively, as follows.

Figure 1: Two-point block method.

First Point

Second Point

The function in these equations will be approximated using Lagrange interpolating polynomial, By following the standard mathematical process, replace the in (4) and (6) with Lagrange interpolation polynomial that interpolates the set of points

Now, introducing the variable substitution and into interpolating polynomial, hence, evaluate these integrals using MAPLE with the limit of the integration from to This will generate the corrector formula for the first point. Similar procedure applied to the in (8) and (10) with Lagrange interpolation polynomial that interpolates the set of points

Next, introduce the variable substitutions and into interpolating polynomial. Again, evaluate these integrals using MAPLE with the limit of the integration from to This will generate the corrector formula for the second point. This two-point one block method is the combination of mode where P is the predictor formula, C is the corrector, and E is the evaluation of the function, The predictor formulae were derived similarly as the corrector formulae but with the order being one less. In this study, we called the proposed predictor-corrector formula as 2PDD4 method. The 2PDD4 formulae were given as follows.

Predictor

Corrector

One-step method will be used at the beginning of the proposed block method in order to get the starting initial points since multistep method needs more than one previous points before generating the remaining values over the interval.

3. Analysis of the Method

In this section, the order, consistency, stability, and convergence of the proposed two-point block method are discussed.

3.1. Order of the Method

The main proposed method of (13)–(16) is classified as a member of the Linear Multistep Method (LMM) which generally can be represented as

The local truncation error (LTE) associated with LMM in (17) is defined in the form of linear difference operator as

Assume that is sufficiently differentiable, so that expanding the terms in (18) using Taylor’s series about the point will give the following expression:

Substituting (19)–(21) into (18) results in

Simplifying this gives

and now we obtained the expression

where

Definition 1. According to Fatunla [7] and Lambert [14], the method in (17) is said to be of order with error constant if

This concept was used to calculate the order and error constant of the proposed formula as stated in (13)–(16). Now, rewrite the corrector formula in matrix difference form as

with

By choosing the calculation gives

Thus, the 2PDD4 satisfied order four with the error constant,

3.2. Consistency of the Method

Definition 2. The linear multistep method is said to be consistent if it has order according to [14].
Since the order of the proposed method is , the method is consistent.

3.3. Stability of the Method

Definition 3. The linear multistep method is zero stable provided that the root of the first characteristics polynomial specified as satisfies and for those roots with , the multiplicity must not exceed two as in Lambert [14] and See et al. [15].
Rewrite (15) to (16) in matrix form as follows:

From (30), the first characteristic polynomial , where

According to Definition 3 and (31), the diagonal two-point one block method is zero stable.

The test equation used in order to obtain the stability polynomial of the two-point block follows the idea from [11] as follows:

The stability polynomial obtained is as follows:

where and

The boundary of the absolute stability region in plane is determined by substituting in the stability polynomial with and for Figure 2 illustrates the regions of the absolute stability for the block method.

Figure 2: Stability region of the two-point diagonal block method.

The stability region shows that the proposed numerical method will be able to produce reasonable results with the given values of the time steps.

3.4. Convergence of the Method

Definition 4. The linear multistep method is convergent if and only if it is consistent and zero stable [14].

Since the consistency and zero stability of the method have been achieved, in conclusion, the proposed two-point block method is convergent.

4. Implementation of the Method

In this study, shooting technique will be adapted throughout the solution process where the calculation starts with deciding a set of initial guesses. The behaviour of shooting method is ‘hit or miss’ the target. Due to that, we attempt to obtain the result as close as possible to the required solution with less number of guessing values. Therefore, choosing the right initial guesses and implementing the best method for refining the guessing value are two important elements in this shooting procedure. Initial values of and need to be guessed for the case of Robin type since none of them are given in (2). Starting with the initial guess, for , then the initial guess of is given explicitly from the first boundary condition as follows:

where and Compute the approximate solutions using the formula in (13) to (16). We obtained the first stopping condition as

where , with If this is sufficiently close to the condition in (35), then the BVPs are solved. Otherwise, update the new set of guessing values and the process continues as described in the following algorithm.

Algorithm of 2PDD4

Step 1. Set TOL and

Step 2. Set and calculate the approximate values, and , using the formula in (13) to (16) with where and depend on the test of convergence at each iteration.

Step 3. If , then repeat Step 2. If , then go to Step 4.

Step 4. If the following stopping condition is fulfilled: , then go to Step 6.
Else, continue to Step 5.

Step 5. Generate the new guessing values, and for based on the previous guess using Newton divided difference interpolation formula. Repeat Step 2.

Step 6. Complete.

The first two estimates were decided to be and based on the consideration in Roberts [16]. The calculation for the corrector part involved the test for convergence. The formulae for the calculation of error and convergence test are defined as follows.

Error

Convergence Test

where r is the number of iterations and is the component of the approximation. result in absolute error test, result in mixed error test, and result in relative error test. All the calculations were done using the C programing code.

5. Numerical Results

In this section, we have applied the algorithm of 2PDD4 to five numerical tested problems to illustrate its accuracy and efficiency. Problem 1 used mixed error test and Problem 2 applied the relative error test whereas Problems 3, 4, and 5 used absolute error test throughout the calculation for obtaining the required result. The following notations are used in the following result.

 MAXE: Maximum absolute error as stated with  AVE: Average absolute error : Step size TOL: Tolerance TS: Total steps at last iteration FCN: Total function call ITN: Total iteration of guess TIME: Time computation in seconds 2PDD4: Direct two-point diagonal block method of order four proposed in this study 2PDAM4: Direct two-step Adams Moulton method of order four as in Phang et al. [17] DAM4: Direct Adams Moulton method of order four as in Majid et al. [18] BP5: Bernoulli polynomials of degree five proposed by Islam and Shirin [2].

Problem 1. One has linear second-order differential equation

with and
Exact solution:
Source: Islam and Shirin [2].

Problem 2. One has nonlinear second-order differential equation

with and
Exact solution:
Source: Duan et al. [3].

Problem 3. One has nonlinear second-order differential equation

with and
Exact solution:
Source: Duan et al. [3].

Problem 4. One has nonlinear second-order differential equation

with and
Exact solution: Source: Bhatta and Sastri [5].

Problem 5. One has nonlinear second-order differential equation

with and
Exact solution:
Source: Lang and Xu [6].
The general equation of Problem 5 given as follows:

arises in applications involving the diffusion of heat generated by positive temperature-dependent sources. According to discussion in Agarwal and O’Regan [19], p.295, if , the diffusion of heat arises in the two different analyses as Joule losses either in electrically conducting solids or in frictional heating. For the first analysis, represents the square of the constant current and the temperature-dependent resistance. Meanwhile, for the second part, represents the square of the constant shear stress and the temperature-dependent fluidity. In addition, if and , then problem in (43) has a unique solution.

Data in Tables 1–4 summarized the result obtained for all tested problems at and with In addition, the values ending with in Tables 1–4 denote the maximum error for those particular results. In order to investigate the competency of this proposed method, 2PDD4 has been compared with 2PDAM4 and DAM4. These two existing methods have the same order as 2PDD4 and will solve the BVP problems directly. The same shooting procedure has been implemented in the algorithm as explained in Section 4.

In Table 1, the maximum error of 2PDD4 is comparable with 2PDAM4 and DAM4 but better than BP5. The total function calls and execution times for 2PDD4 are less compared to 2PDAM4 and DAM4. As can be seen in Table 2, the total function calls and number of initial guessing for 2PDD4 are greater than DAM4 at but comparable in terms of computation time. Besides that, 2PDD4 achieved an excellent accuracy result compared to 2PDAM4. In detail, 2PDD4 has less three- or two-decimal accuracy than 2PDAM4 for solving Problem 2 with and , respectively. Method 2PDD4 needs less total function call and obtained smallest maximum error compared to other methods as displayed in Table 3. Table 4 has shown that 2PDD4 required only one initial guessing as compared to 2PDAM4 and DAM4 for solving Problem 4 at Due to this efficient iterative performance, 2PDD4 gives a superiority result in terms of execution times and total function calls.

Table 1: Comparison of the numerical result for solving Problem 1.

Error at
Error at

2PDD4
2PDAM4
DAM4
2PDD4
2PDAM4
DAM4
BP5

1.5708
4.64e-07
2.60e-07
1.49e-06
2.05e-10
2.05e-10
2.05e-10
0.00e+00

1.7279
6.76e-07
5.81e-07
2.37e-07
1.70e-10
1.71e-10
1.71e-10
2.74e-06

1.8850
8.47e-07
8.32e-07
7.10e-07
1.56e-10
1.63e-10
1.63e-10
8.52e-06

2.0420
6.32e-07
7.43e-07
1.11e-06
1.49e-10
1.66e-10
1.69e-10
6.63e-06

2.1991
5.36e-07
6.73e-07
1.39e-06
1.48e-10
1.78e-10
1.91e-10
4.21e-06

2.3562
3.68e-07
6.25e-07
1.56e-06
1.52e-10
2.00e-10
2.29e-10
1.10e-05

2.5133
3.43e-07
5.78e-07
1.62e-06
1.60e-10
2.31e-10
2.86e-10
5.35e-06

2.6704
1.94e-07
5.28e-07
1.59e-06
1.73e-10
2.72e-10
3.65e-10
6.10e-06

2.8274
2.06e-07
4.72e-07
1.45e-06
1.91e-10
3.25e-10
4.68e-10
9.42e-06

2.9845
6.27e-08
3.96e-07
1.18e-06
2.16e-10
3.92e-10
6.03e-10
2.45e-06

3.1416
9.16e-08
3.05e-07
7.81e-07
2.47e-10
4.76e-10
7.76e-10
0.00e+00

AVE
4.42e-07
5.99e-07
1.31e-06
1.75e-10
2.46e-10
3.16e-10
6.27e-06

TS
6
6
10
51
51
100
-

FCN
96
144
96
114
212
114
-

ITN
3
3
3
1
1
1
-

TIME
0.0030
0.0150
0.0310
0.0150
0.0470
0.0470
-

Table 2: Comparison of the numerical result for solving Problem 2.

Error at
Error at

2PDD4
2PDAM4
DAM4
2PDD4
2PDAM4
DAM4

0.00
4.90e-07
3.05e-05
1.84e-07
0.00e+00
0.00e+00
0.00e+00

0.10
2.61e-07
3.04e-05
4.17e-07
5.32e-12
2.79e-11
8.54e-12

0.20
2.63e-08
3.06e-05
6.61e-07
5.87e-12
1.35e-10
2.23e-11

0.30
2.58e-08
3.13e-05
6.76e-07
6.69e-12
3.24e-10
4.57e-11

0.40
9.67e-09
3.21e-05
6.78e-07
7.76e-12
5.88e-10
7.80e-11

0.50
1.30e-08
3.19e-05
6.68e-07
9.08e-12
9.22e-10
1.19e-10

0.60
7.25e-11
3.15e-05
6.46e-07
1.06e-11
1.32e-09
1.68e-10

0.70
6.15e-09
2.99e-05
6.14e-07
1.24e-11
1.79e-09
2.24e-10

0.80
4.33e-09
2.81e-05
5.73e-07
1.44e-11
2.32e-09
2.88e-10

0.90
3.91e-09
2.47e-05
5.22e-07
1.66e-11
2.91e-09
3.59e-10

1.00
4.63e-09
2.11e-05
4.62e-07
1.90e-11
3.56e-09
4.38e-10

AVE
8.44e-08
3.22e-05
6.10e-07
1.01e-11
1.22e-09
1.55e-10

TS
6
6
10
51
51
100

FCN
120
212
88
228
424
228

ITN
4
4
3
2
2
2

TIME
0.015
0.031
0.015
0.031
0.047
0.046

Table 3: Comparison of the numerical result for solving Problem 3.

Error at
Error at

2PDD4
2PDAM4
DAM4
2PDD4
2PDAM4
DAM4

0.00
0.00e+00
0.00e+00
1.28e-06
0.00e+00
0.00e+00
0.00e+00

0.10
1.68e-08
1.68e-08
1.33e-06
3.72e-13
1.29e-10
3.30e-13

0.20
3.03e-08
3.03e-08
1.38e-06
4.26e-13
7.01e-10
7.56e-14

0.30
4.47e-08
2.36e-08
1.43e-06
5.04e-13
1.64e-09
3.62e-13

0.40
3.12e-08
2.16e-08
1.45e-06
5.99e-13
2.86e-09
9.31e-13

0.50
3.55e-08
2.82e-08
1.46e-06
7.05e-13
4.32e-09
1.59e-12

0.60
2.34e-08
3.73e-08
1.45e-06
8.18e-13
5.96e-09
2.32e-12

0.70
2.23e-08
5.12e-08
1.43e-06
9.36e-13
7.74e-09
3.09e-12

0.80
1.12e-08
6.63e-08
1.40e-06
1.06e-12
9.64e-09
3.90e-12

0.90
7.25e-09
8.41e-08
1.36e-06
1.18e-12
1.16e-08
4.73e-12

1.00
2.96e-09
1.03e-07
1.33e-06
1.31e-12
1.37e-08
5.57e-12

AVE
2.26e-08
4.62e-08
1.53e-06
7.46e-13
5.20e-09
2.05e-12

TS
6
6
10
51
51
100

FCN
32
80
32
114
604
114

ITN
1
1
3
1
1
1

TIME
0.0040
0.0470
0.0470
0.0400
0.0620
0.0780

Table 4: Comparison of the numerical result for solving Problem 4.

Error at
Error at

2PDD4
2PDAM4
DAM4
2PDD4
2PDAM4
DAM4

0.00
0.00e+00
2.07e-06
1.18e-05
0.00e+00
0.00e+00
0.00e+00

0.10
5.30e-07
1.77e-06
1.25e-05
1.27e-11
2.41e-11
9.41e-11

0.20
9.48e-07
1.62e-06
1.36e-05
1.20e-11
8.26e-11
3.66e-10

0.30
1.66e-06
2.03e-06
1.46e-05
1.01e-11
1.71e-10
7.57e-10

0.40
1.33e-06
2.36e-06
1.47e-05
7.54e-12
2.79e-10
1.23e-09

0.50
1.53e-06
2.39e-06
1.44e-05
4.68e-12
4.00e-10
1.75e-09

0.60
1.34e-06
2.39e-06
1.37e-05
1.61e-12
5.33e-10
2.33e-09

0.70
1.38e-06
2.29e-06
1.29e-05
1.62e-12
6.74e-10
2.93e-09

0.80
1.25e-06
2.18e-06
1.20e-05
4.97e-12
8.23e-10
3.58e-09

0.90
1.23e-06
2.03e-06
1.10e-05
8.42e-12
9.80e-10
4.25e-09

1.00
1.14e-06
1.88e-06
1.01e-05
1.20e-11
1.14e-09
4.96e-09

AVE
1.23e-06
2.30e-06
1.41e-05
7.55e-12
4.59e-10
2.00e-09

TS
6
6
10
51
51
100

FCN
32
48
32
114
212
114

ITN
1
4
4
1
1
1

TIME
0.0030
0.0150
0.0150
0.0050
0.0310
0.0460

Overall, a similar number of total steps are observed from all the data generated by 2PDD4 and 2PDAM4 as presented in Tables 1–4. These results justify that both methods have potential to reduce the total step approximately by half when computing the results at two-point simultaneously.

Finally, Table 5 describes the comparison of maximum error obtained in Lang and Xu [6] using the Quintic b-spline collocation method of order four with 2PDD4 method. The first three results demonstrate that the maximum error for both methods is comparable. However, at and , the results showed that 2PDD4 are able to achieve smaller maximum error than the method in [6]. Figures 3, 4, 5, and 6 display clear information about the comparison of the maximum errors versus the total function call as well as time for the numerical results summarized in Tables 1–4. Comparison result for Problem 5 is as depicted in Figure 7.

Table 5: Comparison of MAXE for solving Problem 5.

MAXE
Lang and Xu[6]
MAXE of 2PDD4

5
0.2
5.085000E-03
1.929712E-02

10
0.1
3.405000E-04
5.903279E-04

20
0.05
2.154000E-05
2.478867E-05

40
0.025
1.348000E-06
1.042691E-07

80
0.0125
8.425000E-08
6.295665E-09

(a) Comparison of FCN and MAXE

(b) Comparison of TIME and MAXE

Figure 3: Graph of numerical result for solving Problem 1.

(a) Comparison of FCN and MAXE

(b) Comparison of TIME and MAXE

Figure 4: Graph of numerical result for solving Problem 2.

(a) Comparison of FCN and MAXE

(b) Comparison of TIME and MAXE

Figure 5: Graph of numerical result for solving Problem 3.

(a) Comparison of FCN and MAXE

(b) Comparison of TIME and MAXE

Figure 6: Graph of numerical result for solving Problem 4.

Figure 7: Graph of numerical result for solving Problem 5.

6. Conclusion

This study has presented a diagonal two-point block method formula of order four to solve directly the linear and nonlinear second-order Robin boundary conditions. Overview from the results showed a significant finding that the proposed diagonal block method manages to give faster execution time and less number of total function calls in all tested problems. The proposed method is able to achieve good accuracy when solving the problems. Future work should be devoted to solving third- and fourth-order Robin boundary value problems.

Data Availability

All data generated or analysed from the five tested problems in the article during this study are included in this published article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors gratefully acknowledge the financial support received from Putra Grant (Project Code: GP-IPS/2018/9625100), Universiti Putra Malaysia (UPM), and SLAB Scholarship sponsored by the Ministry of Higher Education (MOHE), Malaysia.

References

T. T. Akano and O. A. Fakinlede, “Numerical computation of Sturm-liouville problem with Robin boundary condition,” in Proceedings of the World Academy of Science, Engineering and Technology, International Journal of Mathematical, Computational, Physical, Electrical and Computer Engineering, vol. 9, pp. 684–688, 2015.

M. S. Islam and A. Shirin, “Numerical solutions of a class of second order boundary value problems on using Bernoulli polynomials,” Applied Mathematics, vol. 2, no. 9, pp. 1059–1067, 2011.

J. Duan, R. Rach, A. Wazwaz, T. Chaolu, and Z. Wang, “A new modified Adomian decomposition method and its multistage form for solving nonlinear boundary value problems with Robin boundary conditions,” Applied Mathematical Modelling, vol. 37, no. 20–21, pp. 8687–8708, 2013.

R. Rach, J.-S. Duan, and A.-M. Wazwaz, “Solution of higher-order, multipoint, nonlinear boundary value problems with high-order robin-type boundary conditions by the adomian decomposition method,” Applied Mathematics & Information Sciences, vol. 10, no. 4, pp. 1231–1242, 2016.

S. K. Bhatta and K. S. Sastri, “Symmetric spline procedures for boundary value problems with mixed boundary conditions,” Journal of Computational and Applied Mathematics, vol. 45, no. 3, pp. 237–250, 1993.

F.-G. Lang and X.-P. Xu, “Quintic B-spline collocation method for second order mixed boundary value problem,” Computer Physics Communications, vol. 183, no. 4, pp. 913–921, 2012.

S. O. Fatunla, “A class of block methods for second order ivps,” International Journal of Computer Mathematics, vol. 55, no. 1-2, pp. 119–133, 1995.

P. S. Phang, Z. A. Majid, and M. Suleiman, “On the solution of two point boundary value problems with two point direct method,” International Journal of Modern Physics: Conference Series, vol. 9, pp. 566–573, 2012.

Z. A. Majid, M. Mughti Hasni, and N. Senu, “Solving second order linear Dirichlet and Neumann boundary value problems by block method,” IAENG International Journal of Applied Mathematics, vol. 43, no. 2, pp. 71–76, 2013.

Z. Omar and O. Adeyeye, “Solving two-point second order boundary value problems using two-step block method with starting and non-starting values,” International Journal of Applied Engineering Research, vol. 11, no. 4, pp. 2407–2410, 2016.

Z. A. Majid, Parallel block methods for solving ordinary differential equations [Ph.D. thesis], Universiti Putra Malaysia, 2004.

I. S. M. Zawawi, Z. B. Ibrahim, F. Ismail, and Z. A. Majid, “Diagonally Implicit Block Backward Differentiation Formulas for Solving Ordinary Differential Equations,” International Journal of Mathematics and Mathematical Sciences, vol. 2012, Article ID 767328, 8 pages, 2012.

N. Zainuddin, Z. B. Ibrahim, and K. I. Othman, “Diagonally implicit block backward differentiation formula for solving linear second order ordinary differential equations,” in Proceedings of the 3rd international conference on fundamental and applied sciences (icfas 2014): innovative research in applied sciences for a sustainable future, pp. 69–75, Kuala Lumpur, Malaysia, 2014.

J. D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley & Sons, New York, NY, USA, 1973.

P. P. See, Z. Abdul Majid, and M. Suleiman, “Three-step block method for solving nonlinear boundary value problems,” Abstract and Applied Analysis, vol. 2014, Article ID 379829, 8 pages, 2014.

J. Roberts, Ordinary differential equations, Prentice-Hall, Inc., Englewood Cliffs, N.J., 1979.

P. S. Phang, Z. A. Majid, and M. Suleiman, “Solving nonlinear two point boundary value problem using two step direct method (menyelesaikan masalah nilai sempadan dua titik tak linear menggunakan kaedah langsung dua langkah,” Journal of Quality Measurement and Analysis, vol. 7, no. 1, pp. 129–140, 2011.

Z. A. Majid, P. P. See, and M. Suleiman, “Solving directly two point non linear boundary value problems using direct Adams Moulton method,” Journal of Mathematics and Statistics, vol. 7, no. 2, pp. 124–128, 2011.

R. P. Agarwal and D. ORegan, An introduction to ordinary differential equations, Springer Science & Business Media, 2008.