International Journal of Analysis and Applications 

Volume 19, Number 6 (2021), 929-948 

URL: https://doi.org/10.28924/2291-8639 

DOI: 10.28924/2291-8639-19-2021-929 

 

 

Received May 15th, 2021; accepted June 24th, 2021; published November 16th, 2021. 

2010 Mathematics Subject Classification. 34A45. 

Key words and phrases. stiff systems; variable step; variable order; Adams family; Milne’s estimate. 

©2021 Authors retain the copyrights 

of their papers, and all open access articles are distributed under the terms of the Creative Commons Attribution License. 

929 

 

A COMPUTATIONAL STRATEGY OF VARIABLE STEP, VARIABLE ORDER FOR 

SOLVING STIFF SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 

 

J. G. OGHONYON*, P. O. OGUNNIYI, I. F. OGBU 

 

Department of Mathematics, Covenant University, Ota, Nigeria 

*Corresponding author: godwin.oghonyon@covenantuniversity.edu.ng  

 

ABSTRACT: This research study focuses on a computational strategy of variable step, variable order (CSVSVO) for 

solving stiff systems of ordinary differential equations. The idea of Newton’s interpolation formula combine with 

divided difference as the basis function approximation will be very useful to design the method. Analysis of the 

performance strategy of variable step, variable order of the method will be justified. Some examples of stiff systems 

of ordinary differential equations will be solved to demonstrate the efficiency and accuracy. 

 

 

NOMENCLATURE 

CSVSVO: errors in CSVSVO for solving test application problem 1, 2 and 3.  

Memployed: approach employed. 

Maxerrors: the magnitude of the maximum errors of CSVSVO.  

ConvCriteria: convergence criteria 

Source of Application Problem I: see [5] for more info 

Source of Application Problem II: see [28] for more info. 

Source of Application Problem III: see [18] for mnore info. 

 

 

https://doi.org/10.28924/2291-8639
https://doi.org/10.28924/2291-8639-19-2021-929


Int. J. Anal. Appl. 19 (6) (2021) 930 

 

1.    INTRODUCTION 

In diverse applied sciences, like chemical kinetics, mass-spring-damper systems, and 

control system analysis, we find systems of differential equations whose analytical solutions 

comprise terms with magnitudes that change at rates that are substantially unlike. For instance, 

whenever the analytical solution includes the terms 𝑒−𝑎𝑡 and 𝑒−𝑏𝑡, with 𝑎,𝑏 > 0, where the 

magnitude of 𝑎 is majorly greater than 𝑏, then 𝑒−𝑎𝑡 decays to zero at extremely quicker rate than 

𝑒−𝑏𝑡 does. In cases of a quickly decaying transient analytical solution, a sure computational 

technique turns unstable except the step length is immoderately small. Explicit techniques 

universally are submitted to this stability control, which necessitates the usage of very small 

step length not only necessarily improve the amount of functions to find the analytic solution, 

and as such stimulates round-off error to spring up, hence, having bounded accuracy and 

efficiency. Implicit techniques, then again, are release of stability limitations and are thus 

favourable for computing stiff systems differential equations [27]. 

The conception of stiff initial value problems can be best valued by studying the 

succeeding general one-dimensional systems with changeless constants: 

𝑦′ = 𝐴𝑦 + 𝐵(𝑥),  𝑦(𝑎) = 𝑦0,        (1) 

where 𝐴 is an 𝑚 × 𝑚 matrix with actual entries and 𝐵(𝑥),𝑦,𝑦′ are 𝑚 − 𝑣𝑒𝑐𝑡𝑜𝑟𝑠. 

The theoretic solutions to (1) is seen as  

 𝑦(𝑥) = ∑ 𝛼𝑖𝑒
𝜆𝑖𝑥𝑐𝑖

𝑚
𝑖=1 + 𝑦𝑝(𝑥),        (2) 

where 𝜆𝑖, 𝑖 = 1(1)𝑚  are the eigenvalues of 𝐴 , with 𝑐𝑖, 𝑖 = 1(1)𝑚  the matching eigenvectors. 

𝑦𝑝(𝑥) is a special solution to (1), and 𝛼𝑖, 𝑖 = 1(1)𝑚 are actual constants that are unambiguously 

determined by the related initial conditions 𝑦(𝑎) = 𝑦0 [14, 18-19]. 

Definition 1:  A solution vector (or solution) of the system (1) on the interval 𝐼 is an 

𝑚 × 1 matrix (or vector) of the form 

𝑦(𝑥) = (

𝑦1(𝑥)

𝑦2(𝑥)
⋮

𝑦𝑚(𝑥)

), 

where the 𝑦𝑖(𝑥) are differentiable functions that gratifies (1) on 𝐼[1] for details. 

Definition 2: Any set {𝛼1}𝑖=1
𝑚 = {(

𝛼1𝑖
𝛼2𝑖
⋮
𝛼𝑚𝑖

)}

𝑖=1

𝑚

of 𝑚 linearly independent solution vectors 

of 𝑦′ = 𝐴𝑦 on an interval 𝐼 is called a fundamental set of solutions of 𝑦′ = 𝐴𝑦 on 𝐼. See [1] for 

more info. 



Int. J. Anal. Appl. 19 (6) (2021) 931 

 

Definition 3: The stiffness ratio S of the system (1) is established as  

  𝑆 = {
|𝜆𝑖|(𝑏−𝑎)

In(TOL)
},         (3) 

where In(TOL) is the exponential logarithm of TOL. 

The stiffness ratio as established by (3) is a standard of the dispersion of the fourth dimension 

constants for (1), and in actual problems, may be of the order of 108. See [8] for more info. 

Definition 4: The initial value problem (1) is stated to be stiff if it gratifies (4) and (5); i.e., 

whenever 

(i) 𝑟𝑒(𝜆𝑖) < 0,𝑖 = 1(1)𝑚, and  

(ii) the stiffness ratio 𝑠 > 1 .  

Nevertheless, it should be observed that this is a quite general resolution with respect to 

mathematics. Stiffness takes place whenever the step length is limited by stability, rather than 

order, conditions. See [8] for more info. 

Definition 5:  The initial value problem  

𝑦′ = 𝑓(𝑥,𝑦),  𝑦(𝑎) = 𝑦0,  𝑦 = (𝑦1,𝑦2,… ,𝑦𝑚)
𝑇,𝑦0 = (𝜂1,𝜂2,… ,𝜂𝑚)

𝑇 

is stated to stiff oscillatory whenever the eigenvalues 𝜆𝑖 = 𝑢𝑖 + 𝑗𝑣𝑖, 𝑖 = 1(1)𝑚 of the Jacobian 𝐽 =

(
𝜕𝑓

𝜕𝑑𝑦
) have the succeeding attributes: 

𝑢𝑖 < 0,, 𝑖 = 1(1)𝑚,     

 Max
1≤𝑖≤𝑚

|𝑢𝑖| > min
1≤𝑖≤𝑚

|𝑢𝑖|,        

or whenever the stiffness ratio gratifies  

max
1≤𝑖≤𝑚

|
𝑢𝑖
𝑢𝑗
| > 1 

and   

|𝑢𝑖| < |𝑣𝑖| 

for at least single pair of 𝑖 ∈ 1 ≤ 𝑖 ≤ 𝑚. See [8] for more info. 

Definition 6: Stiffness occurs when stability requirements, rather than those of accuracy 

constrain the step length. See [18-19] for details. 

Definition 7: Stiffness occurs when some components of the solution decay much more 

rapidly than others. See [18-19] for details. 

Authors have contributed immensely towards solving stiff systems of ordinary 

differential equations using diverse strategies. [9] implemented the extensions of the predictor-

corrector method for the solution of systems of ordinary differential equations. [11] formulated 



Int. J. Anal. Appl. 19 (6) (2021) 932 

 

the resultant of variable mesh size on the constancy of multi-step methods. [12] constructed the 

constancy and convergency of variable order multi-step methods. [14] developed the diagonally 

implicit block backward differentiation formula with optimal stability properties for stiff 

ordinary differential equations. [15] launched a varying-step, varying-order multistep method 

for the numerical solution of ordinary differential equations. [16] designed the algorithms for 

changing the  step size. [17] worked on changing step-size in the integration of differential 

equations using modified divided differences. [21-25] developed and implemented a variable 

step, variable step size with same order for solving ordinary differential equations. [26] derived 

the constancy, consistence and convergency of varying K-step methods for numerical 

integration of large-systems of ordinary differential equations.  This research is designed to 

extend the idea of [21-25] by inventing variable step, variable order together with variable step 

size for solving stiff systems of ordinary differential equations. [18-19] specifies that this idea 

leads to better efficiency and accuracy and as well bypass theorem 4.  

Nevertheless, nonstiff algorithmic program possess a limited region of absolute stability 

(RAS), whilst stiff algorithmic program possess unlimited RAS. This describes why stiff 

algorithmic program adapt the usage of large step length beyond the transient (nonstiff) phase. 

In the transient phase (boundary layer), automatic codification seeks to name the optimal mesh 

size that holds the local truncation error inside the bound of the observed accuracy. Again, for 

the transient phase, accuracy necessities oblige a numeric integrator to accept a mesh size of the 

order of the littlest time constant. Beyond the transient phase (i.e., in the stiff phase) the increase 

of disseminated errors (instability) checks the selection of mesh size whenever a nonstiff 

technique is followed, since stability limitations are autonomous of the accuracy necessities. For 

the transient phase, there is ever the need to accept fairly large mesh size, but this is frequently 

bounded by stability conditions. See [8] for more info. 

Theorem 1: The set {𝛼1}𝑖=1
𝑚 = {(

𝛼1𝑖
𝛼2𝑖
⋮
𝛼𝑚𝑖

)}

𝑖=1

𝑚

is linearly independent if and only if the 

Wronskian 

𝑊({𝛼1}𝑖=1
𝑚 ) = |𝛼1𝛼2 …𝛼𝑚| = |

𝛼11𝛼12 …𝛼1𝑚
𝛼21𝛼22 …𝛼2𝑚
⋮           ⋮ ⋯ ⋮

𝛼𝑚1𝛼𝑚2 …𝛼𝑚𝑚

| ≠ 0  

See [1] for more info. 



Int. J. Anal. Appl. 19 (6) (2021) 933 

 

Theorem 2:  Let 𝑆 = {𝛼1}𝑖=1
𝑚 = {(

𝛼1𝑖
𝛼2𝑖
⋮
𝛼𝑚𝑖

)}

𝑖=1

𝑚

be a set of 𝑚  linearly independent 

solutions of 𝑦′ = 𝐴𝑦. Then every solution of 𝑦′ = 𝐴𝑦 is a linear combination of these solutions. 

See [1] for more info. 

Theorem 3 (Existence and Uniqueness): Assume that  

each of the functions 

𝑓1(𝑥,𝑦1,𝑦2,…,𝑦𝑚), 𝑓2(𝑥,𝑦1,𝑦2,…,𝑦𝑚),…,𝑓𝑚(𝑥,𝑦1,𝑦2,…,𝑦𝑚) 

and the partial derivatives 
𝜕𝑓1

𝜕𝑥1
,
𝜕𝑓2

𝜕𝑥2
,…,

𝜕𝑓𝑚

𝜕𝑥𝑚
 are continuous in a region 𝑅 containing the point 

(𝑥,𝑦1,𝑦2,…,𝑦𝑚). Then, the initial-value problem 

{
 
 

 
 

𝑦1
′ = 𝑓1(𝑥,𝑦1,𝑦2,…,𝑦𝑚)

𝑦2
′ = 𝑓2(𝑥,𝑦1,𝑦2,…,𝑦𝑚) 

⋮
𝑦𝑚
′ = 𝑓𝑚(𝑥,𝑦1,𝑦2,…,𝑦𝑚)

𝑦1(𝑥0) = 𝑡1,… ,𝑦𝑚(𝑥0) = 𝑡𝑚

        (4) 

has a unique solution 

  {

𝑦1 = ∅1(𝑥)

𝑦2 = ∅2(𝑥)
⋮

𝑦𝑚 = ∅𝑚(𝑥)

                    (5) 

on the interval 𝐼 contanining 𝑥 = 𝑥0. See [1, 3] for more info. 

  Theorem 4 (Dahlquist Barrier Theorem): An A − stable linear multistep method  

• must be implicit and  

• The most accurate A-stable linear multistep method is the trapezoidal scheme of order 

𝑝 = 2 and error constant 𝑐3 = −
1

12
 . See [8] for more info. 

The Dahlquist Barrier Theorem 4 can be outwitted by accepting unconventional numeric 

integrators, some of which are 

• nonlinear multistep schemes, 

• multiderivative multistep schemes, 

• exponentially fitting, and  

• extrapolation process.  

The variable step, variable order of the predictor-corrector algorithmic program came forth 

as result of the broad computational experience that holds throughout the years. This VSVO-P-



Int. J. Anal. Appl. 19 (6) (2021) 934 

 

CAP is fundamental to high level efficiency and accuracy with the potential to change 

automatically not exclusively the step length but as well the order (and thus the step number) of 

the techniques utilized. Algorithmic program with such a potential are recognized as variable 

step, variable order or VSVO, algorithmic program. Apart from been built to handle nonstiff 

initial value problems, though various subsisting VSVO codification admit alternative for stiff 

systems. The necessary elements of VSVO algorithmic programs are: 

• a family of methods, 

• a starting procedure, 

• a local error estimator, 

• a strategy for determining when to vary step length and/ or order, and  

• a technique for varying step length and/ or order 

• a written softcode in any mathematical packages is required if manual computation is 

very tedious. 

• a special basis function approximation for estimating stiff systems is necessary if the 

desired accuracy and efficiency is not achieved.  

In addition, the convergence attributes of predictor-corrector methods antecedently proved 

on the presumption of changeless step length and changeless order still maintain in a VSVO 

algorithmic program conceptualizations. The results of [11-12] indicate that a VSVO algorithmic 

program established on Adams-Bashforth-Moulton pair in 𝑃𝐸𝐶𝐸  mode (ABM) with step-

varying attained by a variable coefficient technique is ever convergent  ( as the maximal step 

length used in the time interval of integration inclines to zero). Whenever an interpolatory 

technique is utilized then convergence is ensured whenever the step/order-varying technique is 

such that there subsists a changeless 𝑁 such that in any 𝑁 sequential steps there are ever 𝑘 𝑠𝑡𝑒𝑝𝑠 

of changeless length considered by the same 𝑘𝑡ℎ − 𝑜𝑟𝑑𝑒𝑟 ABM method, for some value of 𝑘. 

These results stress so far once again that variable coefficient technique, although in general is 

more costly and cumbersome to carry out, are essentially more effective than interpolatory 

techniques. See [18-19] for more details. 

The succeeding sections will demonstrate the usefulness of these strategies.  

 

 



Int. J. Anal. Appl. 19 (6) (2021) 935 

 

II. MATERIALS AND METHODS 

[4] was evidently the first author to propose a “placid” conversion of a step-size ℎ to 

novel step-size 𝑤ℎ. [9, 15] unfolded his ideas:  we study an arbitrary grid (𝑥𝑛) and announce the 

step sizes by ℎ𝑛 = 𝑥𝑛+1 − 𝑥𝑛. We presume that estimations 𝑦𝑖 to 𝑦(𝑥𝑖) are recognized for 𝑖 = 𝑛 −

𝑘 + 1,…,𝑛 and insert 𝑓𝑖(𝑥𝑖,𝑦𝑖) and announce 𝑝(𝑥) as the multinomial which interpolates the 

values (𝑥𝑖,𝑓𝑖) for 𝑖 = 𝑛 − 𝑘 + 1,…,𝑛. Utilizing Newton’s interpolation formula we get 

𝑝(𝑥) = ∑ ∏ (𝑥 − 𝑥𝑛−𝑖)𝛿
𝑖𝑓[𝑥𝑛,𝑥𝑛−1,…,𝑥𝑛−𝑖]

𝑖−1
𝑖=0

𝑘−1
𝑖=0     (6) 

where the divided differences 𝛿𝑖𝑓[𝑥𝑛,𝑥𝑛−1,…,𝑥𝑛−𝑖] are 

specified algorithmic by  

𝛿𝑖𝑓[𝑥𝑛] = 𝑓𝑛 

𝛿𝑖𝑓[𝑥𝑛,𝑥𝑛−1,…,𝑥𝑛−𝑖] =
𝛿𝑖−1𝑓[𝑥𝑛,𝑥𝑛−1,…,𝑥𝑛−𝑖+1]−𝛿

𝑖−1𝑓[𝑥𝑛−1,…,𝑥𝑛−𝑖]

𝑥𝑛−𝑥𝑛−𝑖
                                  (7) 

It is virtual to rescript (6) as  

𝑝(𝑥) = ∑ ∏
𝑥−𝑥𝑛−𝑖

𝑥𝑛+1−𝑥𝑛−𝑖
∙ Φ𝑖

∗(𝑛)𝑖−1𝑖=0
𝑘−1
𝑖=0 ,    (8) 

where  

Φ𝑖
∗(𝑛) = ∏ (𝑥𝑛+1 − 𝑥𝑛−𝑖) ∙ 𝛿

𝑖𝑓[𝑥𝑛,𝑥𝑛−1,…,𝑥𝑛−𝑖]
𝑖−1
𝑖=0       (9) 

We immediately specify the estimations to 𝑦(𝑥𝑛+1) by 

𝑦𝑛+1 = 𝑦𝑛 + ∫ 𝑝(𝑥)𝑑𝑥
𝑥𝑛+1
𝑥𝑛

     (10) 

Replacing equation (6) into (10) we have 

𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑛 ∑ 𝑔𝑖(𝑛)Φ𝑖
∗(𝑛)𝑘−1𝑖=0     (11) 

with  

𝑔𝑖(𝑛) =
1

ℎ𝑛
∫ ∏

𝑥−𝑥𝑛−𝑖

𝑥𝑛+1−𝑥𝑛−𝑖
𝑑𝑥𝑖−1𝑖=0

𝑥𝑛+1
𝑥𝑛

.    (12) 

Equation (11) is the elongation of the explicit Adams method 

(1) to variable step sizes [13]. 

The variable step size implicit Adams methods can be inferred likewise. We assume 

𝑝∗(𝑥) be the multinomial of degree 𝑘 that interpolates (𝑥𝑖,𝑓𝑖) for 𝑖 = 𝑛 − 𝑘 + 1 (the value 𝑓𝑛+1 =

𝑓(𝑥𝑛+1,𝑦𝑛+1) comprises the unknown physical solution 𝑦𝑛+1). Once more, employing Newton’s 

interpolation 

formula, we get 

𝑝∗(𝑥) = 𝑝(𝑥) + ∏ (𝑥 − 𝑥𝑛−𝑖) ∙ 𝛿
𝑘𝑓[𝑥𝑛+1,𝑥𝑛,… ,𝑥𝑛−𝑘+1]

𝑖−1
𝑖=0 . 

The numeric solution specified by  



Int. J. Anal. Appl. 19 (6) (2021) 936 

 

𝑦𝑛+1 = 𝑦𝑛 + ∫ 𝑝
∗(𝑥)𝑑𝑥

𝑥𝑛+1
𝑥𝑛

, 

is established immediately by 

𝑦𝑛+1 = 𝑝𝑛+1 + ℎ𝑛𝑔𝑘(𝑛)Φ𝑘(𝑛 + 1)              (13) 

where 𝑝𝑛+1 is the numeric estimation got by the explicit Adams method 

𝑝𝑛+1 = 𝑦𝑛 + ℎ𝑛 ∑ 𝑔𝑖(𝑛)Φ𝑖
∗(𝑛)

𝑘−1

𝑖=0

 

and where  

Φ𝑘(𝑛 + 1) = ∏ (𝑥𝑛+1 − 𝑥𝑛−𝑖) ∙ 𝛿
𝑘𝑓[𝑥𝑛+1,𝑥𝑛,… ,𝑥𝑛−𝑘+1]

𝑖−1
𝑖=0 .      (14) 

Defining the reoccurrence relations for 𝑔𝑖(𝑛),Φ𝑖(𝑛) and Φ𝑖
∗(𝑛). 

The price of calculating consolidation coefficients is the greatest weakness to allowing arbitrary 

fluctuations in the step Size [16]. 

The values Φ𝑖
∗(𝑛)(𝑖 = 0,…,𝑘 − 1)  and Φ𝑘(𝑛 + 1)  can be calculated expeditiously with 

reoccurrence relations 

Φ0(𝑛) = Φ𝑖
∗(𝑛) = f𝑛 

Φ𝑖+1(𝑛) = Φ𝑖(𝑛)− Φ𝑖
∗(𝑛 − 1)              (15) 

Φ𝑖
∗(𝑛) = 𝛽𝑖(𝑛)Φ𝑖(𝑛), 

which are an instant effect of (9) and (14). The coefficients 

𝛽𝑖(𝑛) = ∏
𝑥𝑛+1 − 𝑥𝑛−𝑖
𝑥𝑛 − 𝑥𝑛−𝑖−1

𝑑𝑥
𝑖−1

𝑖=0
 

can be computed by  

𝛽0(𝑛) = 1,  𝛽𝑖(𝑛) = 𝛽𝑖−1(𝑛)
𝑥𝑛+1−𝑥𝑛−𝑖+1

𝑥𝑛−𝑥𝑛−1
. 

The computation of the coefficients 𝑔𝑖(𝑛) is wilier [17].  We  

bring in the 𝑞 − 𝑓𝑜𝑙𝑑 integral  

𝑐𝑖𝑞(𝑥) =
(𝑞−1)!

ℎ𝑛
𝑞 ∫ …

𝑥

𝑥𝑛
∫ ∏

𝜖0−𝑥𝑛−𝑖

𝑥𝑛+1−𝑥𝑛−𝑖
𝑑𝜖0 …

𝑖−1
𝑖=0

𝜖𝑞−1
𝑥𝑛

𝑑𝜖𝑞−1           (16) 

and remark that  

𝑔𝑖(𝑛) = 𝑐𝑖𝑞(𝑥𝑛+1) [13]. 

 

A. Theoretical Analysis of the Method 

Definition 8: The order of the operator 𝐿ℎ is the highest 𝑟 such that whenever 𝑦(𝑥) possesses a 

continuous 

(𝑟 + 1)𝑡ℎ the derivative, then  



Int. J. Anal. Appl. 19 (6) (2021) 937 

 

𝐿ℎ(𝑦(𝑥)) = 0(ℎ
𝑟+1).      (17) 

Whenever we presume a continuous (𝑟 + 2)𝑡ℎ derivative for 𝑦, then we can replace the Taylor’s 

series for 𝑦 and 𝑦′ with 0(ℎ𝑟+1) remains. Whenever the terms in ℎ0,ℎ2,…,ℎ𝑟+1 are  

assembled unitedly, we will arrive at  

𝐿ℎ(𝑦(𝑥)) = ∑ 𝐶𝑞
𝑟+1

𝑞=0
ℎ𝑞𝑦(𝑞)(𝑥) + 0(ℎ𝑟+1) 

where  

𝐶𝑞 = {
∑ 𝛼𝑖
𝑘
𝑖=0                                                 𝑞 = 0

∑ [
(−𝑖)𝑞

𝑞!
𝛼𝑖 −

(−𝑖)𝑞−1

(𝑞−1)!
𝛽𝑖]        𝑞 > 0

𝑘
𝑖=0

     (18) 

The linear equations 𝐶𝑞 = 0,𝑞 ≤ 𝑟, are the equations which decides an 𝑟𝑡ℎ order method.  

Given the number of truncation error put in for each one step is  

𝐶𝑟+1

∑ 𝛽𝑖
𝑘
𝑖=0

ℎ𝑟+1𝑦(𝑟+1) + 0(ℎ𝑟+2). 

Therefore the natural standization is to take 

∑ 𝛽𝑖
𝑘
𝑖=0 = 1 [10].       (19) 

Theorem 5: Whenever the multi-step method (29) is unchanging and of order 1  so it is 

convergent. Whenever the method (29) is unchanging and of order 𝑝 so it is convergent of order 

𝑝 [13]. 

Theorem 6:  The multi-step method (29) is of order 𝑝, whenever one and only of the 

following tantamount precondition is true conditions is met: 

(i) ∑ 𝛼𝑖
𝑘
𝑖=0 = 0 and ∑ 𝛼𝑖

𝑘
𝑖=0 𝑖

𝑞 = 𝑞∑ 𝛽𝑖𝑖
𝑞−1𝑘

𝑖=0  for 𝑞 = 1,…,𝑝; 

(ii) 𝜚(𝑒ℎ)− ℎ𝜎(𝑒ℎ) = 𝑂(ℎ𝑝+1) for ℎ → 0; 

(iii) 
𝜚(𝜍)

𝑙𝑜𝑔𝜍
− 𝜎(𝜍) = 𝑂((𝜍 − 1)𝑝) for 𝜍 → 1. 

Where the linear difference operator 𝐿 specified by 

𝐿(𝑦,𝑥,ℎ) = ∑ (𝛼𝑖𝑦(𝑥 + 𝑖ℎ) − h𝛽𝑖𝑦
′(𝑥 + 𝑖ℎ))𝑘𝑖=0  [13].     (20) 

Proof 

Enlarging 𝑦(𝑥 + 𝑖ℎ) and 𝑦′(𝑥 + 𝑖ℎ) using Taylor series and inserting the truncated series. 

Put the Taylor series expansion into (20) gives 

𝐿(𝑦,𝑥,ℎ) = ∑ (∑
𝑖𝑞

𝑞!𝑞≥0
ℎ𝑞𝑦(𝑞)(𝑥) − ℎ𝛽𝑖 ∑

𝑖𝑟

𝑟!𝑟≥0
ℎ𝑟𝑦(𝑟+1)(𝑥))𝑘𝑖=0    (21) 

      𝑦(𝑥)∑ 𝛼𝑖
𝑘
𝑖=0 + ∑

ℎ𝑞

𝑞!𝑞≥1
𝑦(𝑞)(𝑥)(∑ 𝛼𝑖𝑖

𝑞 − 𝑞∑ 𝛽𝑖𝑖
𝑞−1𝑘

𝑖=0
𝑘
𝑖=0 ). 



Int. J. Anal. Appl. 19 (6) (2021) 938 

 

This means the par of precondition (i) with 𝐿(𝑦,𝑥,ℎ) = 𝑂(ℎ𝑝+1) for entirely sufficient degree of 

regular functions 𝑦(𝑥). 

We will continue to express that the three precondition of proposition 1 are tantamount. The 

individuality operator  

𝐿(𝑒𝑥𝑝,0,ℎ) = 𝜚(𝑒ℎ) − ℎ𝜎(𝑒ℎ) 

where 𝑒𝑥𝑝 announces the exponential function, unitedly with  

𝐿(𝑒𝑥𝑝,0,ℎ) = ∑ 𝛼𝑖
𝑘
𝑖=0 + ∑

ℎ𝑞

𝑞!𝑞≥1
(∑ 𝛼𝑖𝑖

𝑞 − 𝑞∑ 𝛽𝑖𝑖
𝑞−1𝑘

𝑖=0
𝑘
𝑖=0 ), 

which succeeds from (21), proves the par of the preconditions (i) and (ii). 

By use of the translation 𝜍 = 𝑒ℎ (𝑜𝑟 ℎ = 𝑙𝑜𝑔𝜍) precondition (ii) can be spelt in the pattern 

𝜚(𝜍) − 𝑙𝑜𝑔𝜍.𝜎(𝜍) = 𝑂((𝑙𝑜𝑔𝜍)𝑝+1) for 𝜍 → 1. 

But this precondition is par to (iii), since  

𝑙𝑜𝑔𝜍 = (𝜍 − 1) + 𝑂((𝜍 − 1)2) for 𝜍 → 1. See [13] for more info. 

 

B. Convergence 

Convergence for variable step size Adams method was first considered by [26]. In order 

to show convergence for the general case we present the vector 𝑌𝑛 = (𝑦𝑛+𝑘−1,… ,𝑦𝑛+1,𝑦𝑛)
𝑇. The 

method  

𝑦𝑛+𝑘 + ∑ 𝛼𝑖𝑦𝑛+𝑖 = ℎ𝑛+𝑘−1
𝑘−1
𝑖=0 ∑ 𝛽𝑖𝑛𝑓𝑛+𝑖

𝑘
𝑖=0     (22) 

then turns tantamount to  

𝑌𝑛+1 = (𝐴𝑛 ⊗ 𝐼)𝑌𝑛 + ℎ𝑛+𝑘−1Φ𝑛(𝑥𝑛,𝑌𝑛,ℎ𝑛) ,    (23) 

where 𝐴𝑛 is established by  























=

−−− −

01

001

00...01

......
,0,1,1

A 

 nnnk

n
       (24) 

the comrade matrix and 

Φ𝑛(𝑥𝑛,𝑌𝑛,ℎ𝑛) = (𝑒1 ⊗ 𝐼)Ψ𝑛(𝑥𝑛,𝑌𝑛,ℎ𝑛). 

The value Ψ = Ψ𝑛(𝑥𝑛,𝑌𝑛,ℎ𝑛) is specified without any doubt by 

Ψ = 𝑦𝑛+𝑘 + ∑ 𝛽𝑖𝑛𝑓(𝑥𝑛+𝑖,𝑦𝑛+𝑖) + 𝛽𝑘𝑛𝑓(𝑥𝑛+𝑘,ℎΨ −
𝑘−1
𝑖=0 ∑ 𝛼𝑖𝑛𝑦𝑛+𝑖

𝑘
𝑖=0 . 

We advance by announcing  



Int. J. Anal. Appl. 19 (6) (2021) 939 

 

𝑌(𝑥𝑛) = (𝑦(𝑥𝑛+𝑘−1),…,𝑦(𝑥𝑛+1),𝑦(𝑥𝑛))
𝑇 

the precise values to be estimated by 𝑌𝑛 . The convergence theorem can immediately be 

developed as succeeds. See [13] for more info. 

Theorem 7: Assume that  

• the method (27) is stable of order 𝑝, and has bounded coefficients 𝛼𝑖𝑛 and 𝛽𝑖𝑛; 

• the starting values satisfy ‖𝑌(𝑥𝑛) − 𝑌0‖ = 𝑂(ℎ0
𝑝
); 

• the step size ratios are bounded (
ℎ𝑛

ℎ𝑛−1
≤ Ω). 

Then the method is convergent of order 𝑝, i.e., for each differential equation  

𝑦′ = 𝑓(𝑥,𝑦),𝑦(𝑥0) = 𝑦0 

with 𝑓  sufficiently differentiable the global error satisfies 

‖𝑌(𝑥𝑛) − 𝑦𝑛‖ ≤ 𝐶ℎ
𝑝   for 𝑥𝑛 ≤ 𝑥, 

where ℎ = 𝑚𝑎𝑥ℎ𝑖. See [13] for more info. 

Proof 

Because the approach has order p and the physical coefficients and step-size ratios are 

restricted, the expression becomes 

𝑦(𝑥𝑛+𝑘)+ ∑ 𝛼𝑖𝑛𝑦(𝑥𝑛+𝑖) − ℎ𝑛+𝑘−1
𝑘−1
𝑖=0 ∑ 𝛽𝑖𝑛𝑦

′(𝑥𝑛+𝑖) =
𝑘
𝑖=0 𝑂(ℎ𝑛

𝑝+1
), 

We justify that the truncated-local error  

𝛿𝑛+1 = 𝑌(𝑥𝑛+1)− (𝐴𝑛 ⊗ 𝐼)𝑌(𝑥𝑛) − ℎ𝑛+𝑘−1Φ𝑛(𝑥𝑛,𝑌(𝑥𝑛),ℎ𝑛)   (25) 

gratifies   

𝛿𝑛+1 = 𝑂(ℎ𝑛
𝑝+1

).        (26) 

Deducting (23) from (25) we get  

𝑌(𝑥𝑛+1)− 𝑌𝑛+1 = (𝐴𝑛 ⊗ 𝐼)𝑌(𝑥𝑛) − 𝑌𝑛) +ℎ𝑛+𝑘−1(Φ𝑛(𝑥𝑛,𝑌(𝑥𝑛),ℎ𝑛)− Φ𝑛(𝑥𝑛,𝑌𝑛,ℎ𝑛)) + 𝛿𝑛+1 

and by induction it succeeds that  

𝑌(𝑥𝑛+1) − 𝑌𝑛+1

= (𝐴𝑛 ….𝐴0) ⊗ 𝐼)𝑌(𝑥0) − 𝑌0) + ∑ ℎ𝑖+𝑘−1
𝑛

𝑖=0
(𝐴𝑛 ….𝐴𝑖+1)⊗ 𝐼)(Φ𝑖(𝑥𝑖,𝑌(𝑥𝑖),ℎ𝑖)

− Φ𝑖(𝑥𝑖,𝑌𝑖,ℎ𝑖)) + ∑ (𝐴𝑛 ….𝐴𝑖+1) ⊗ 𝐼)
𝑛

𝑖=0
𝛿𝑖+1. 

As in the proof of theorem 1, we derive that the Φ𝑛 gratifies a uniform Lipschitz precondition 

with respect to 𝑌𝑛. This unitedly with stability and (26), means that  

‖𝑌(𝑥𝑛+1) − 𝑌𝑛+1‖ ≤ ∑ ℎ𝑖+𝑘−1
𝑛

𝑖=0
𝐿‖𝑌(𝑥𝑖) − 𝑌𝑖‖ + 𝐶1ℎ

𝑝. 



Int. J. Anal. Appl. 19 (6) (2021) 940 

 

To figure out this difference, we bring in the succession { 𝑛} specified as  

0 = ‖𝑌(𝑥0) − 𝑌0‖,   

𝑛+1 = ∑ ℎ𝑖+𝑘−1
𝑛
𝑖=0 𝐿 𝑖 + 𝐶1ℎ

𝑝.    (27) 

A simple induction statement proves that  

‖𝑌(𝑥𝑛) − 𝑌𝑛‖ ≤ 𝑛.                  (28) 

From (27) we get for 𝑛 ≥ 1 

𝑛+1 = 𝑛 + ℎ𝑛+𝑘−1𝐿 𝑛 ≤ exp (𝐿ℎ𝑛+𝑘−1) 𝑛 

so that in addition 

𝑛 ≤ exp(𝑥 − 𝑥0)𝐿) 1 = exp(𝑥 − 𝑥0)𝐿) ∙ (ℎ𝑘−1𝐿‖𝑌(𝑥0) − 𝑌0‖ + 𝐶1ℎ
𝑝). 

The inequality unitedly with (28) completes the proof of theorem 7. See [13] for more info. 

 C. Implementing the Convergence Criteria of Variable Step, Variable Order 

The use of Milne’s estimate for the principal local truncation error necessitate that 

predictor-corrector method possess like order. This is attained by accepting the predictor to be a 

𝑘 − 𝑠𝑡𝑒𝑝 Adams Bashforth method and the corrector to be a (𝑘 − 1) − 𝑠𝑡𝑒𝑝 Adams-Moulton, 

both then possess order 𝑝 = 𝑘.  The 𝑘 − 𝑠𝑡𝑒𝑝 𝑘𝑡ℎ order ABM pair is therefore 

 𝑦𝑛+1 = 𝑦𝑛 + ℎ∑ 𝛾𝑖
∗∇𝑖𝑓𝑛,

𝑘−1
𝑖=0      𝑝

∗ = 𝑘,          𝐶𝑘+1
∗ = 𝛾𝑘

∗ 

       𝑦𝑛+1 = 𝑦𝑛 + ℎ∑ 𝛾𝑖∇
𝑖𝑓𝑛+1,

𝑘−1
𝑖=0     𝑝 = 𝑘,           𝐶𝑘+1 = 𝛾𝑘    

}𝑘 = 1,2,…    (29) 

Whenever we imagine (29) being employed in 𝑃(𝐸𝐶)𝜇𝐸1−𝑡 mode then, in the second of (29), 

 𝑦𝑛+1 will be substituted by 𝑦𝑛+1
[𝑣+1]

, and the one value  𝑓𝑛+1 on the right side by 𝑓𝑛+1
[𝑣]

, the leftover 

values  𝑓𝑛−𝑗  being substituted by 𝑓𝑛−𝑗
[𝜇−𝑡]

, 𝑗 = 0,1,…,𝑘 − 1. We can surmount this trouble by 

defining ∇𝑣
𝑖 𝑓𝑛+1

[𝜇]
 to be ∇𝑖𝑓𝑛+1

[𝜇]
 with the one value 𝑓𝑛+1

[𝜇]
 substituted by 𝑓𝑛+1

[𝑣]
. That is, 

∇𝑣
𝑖 𝑓𝑛+1

[𝜇]
= ∇𝑖𝑓𝑛+1

[𝜇]
+ 𝑓𝑛+1

[𝑣]
− 𝑓𝑛+1

[𝜇]
       (30) 

We may rewrite (30) in the form 

∑ (𝛾𝑖∇
𝑖𝑓𝑛+1
[𝜇]

− 𝛾𝑖
∗∇𝑖𝑓𝑛

[𝜇]
) = 𝛾𝑘−1

∗ 𝑓𝑛+1
[𝜇]𝑘−1

𝑖=0          (31) 

where the notational system is specified by (29). 

We immediately employ the pair (29) in 𝑃(𝐸𝐶)𝜇𝐸1−𝑡 mode, and utilize the structure of the 

Adams methods to formulate a form of ABM method which is computationally commodious 

and frugal. The mode is specified by  

𝑃:  𝑦𝑛+1
[0]

= 𝑦𝑛
[𝜇]
+ ℎ∑ 𝛾𝑖

∗∇𝑖𝑓𝑛
[𝜇−1]𝑘−1

𝑖=0        (32) 



Int. J. Anal. Appl. 19 (6) (2021) 941 

 

(𝐸𝐶)𝜇   
𝑓𝑛+1
[𝑣]

= 𝑓(𝑥𝑛+1,𝑦𝑛+1
[𝑣]
)

𝑦𝑛+1
[𝑣+1]

= 𝑦𝑛
[𝜇]
+ ℎ∑ 𝛾𝑖∇𝑣

𝑖 𝑓𝑛+1
[𝜇−1]𝑘−1

𝑖=0

}   𝑣 = 0,1,…,𝜇 − 1     (33) 

        (𝐸1−𝑡   𝑓𝑛+1
[𝜇]

= 𝑓(𝑥𝑛+1,𝑦𝑛+1
[𝜇]
)  

whenever 𝑡 = 0. 

To employ the Milne estimate, we require calculating 𝑦𝑛+1
[𝜇]

− 𝑦𝑛+1
[0]

. Deducting (32) from (33) 

with  𝑣 = 𝜇 − 1 establishes 

𝑦𝑛+1
[𝜇]

− 𝑦𝑛+1
[0]

= ℎ∑ (𝛾𝑖∇𝜇−1
𝑖 𝑓𝑛+1

[𝜇−𝑡]
− 𝛾𝑖

∗∇𝑖𝑓𝑛
[𝜇−𝑡]

)
𝑘−1

𝑖=0
 

= ℎ𝛾𝑖
∗∇𝜇−1

𝑘 𝑓𝑛+1
[𝜇−𝑡]

. 

Since 𝐶𝑘+1
∗ = 𝛾𝑘

∗   and  𝐶𝑘+1 = 𝛾𝑘, the Milne estimate  

𝑊 =
          𝐶𝑝+1

        𝐶𝑝+1
∗ −  𝐶𝑝+1

 

for the principal local truncation error at 𝑥𝑛+1 (which we shall denote by 𝑇𝑛+1) is established by  

𝑇𝑛+1 =
          𝐶𝑝+1

        𝐶𝑝+1
∗ − 𝐶𝑝+1

( 𝑦𝑛+1
[𝜇]

− 𝑦𝑛+1
[0]
) =

𝛾𝑘

𝛾𝑘
∗−𝛾𝑘

= ℎ𝛾𝑖
∗∇𝜇−1
𝑘 𝑓𝑛+1

[𝜇−𝑡]
. 

Whenever 𝛾𝑘
∗ − 𝛾𝑘 = 𝛾𝑘−1

∗ , wherefrom  

𝑇𝑛+1 == ℎ𝛾𝑖
∗∇𝜇−1
𝑘 𝑓𝑛+1

[𝜇−𝑡]
. See [2, 6-7, 18-19, 21-25] for more info. 

 

III. Practical Examples of Stiff Systems of First Order ODEs 

We are interested with the computation interpretation of these attributes. 

A problem is stiff whenever the analytical solution being looked for changes tardily, but there 

are close solutions that change speedily, so the numerical approach must consider small steps 

size to get acceptable results. See [20] for more info. 

Application Problem 1 

An Engineering Example: In chemical engineering, a complicated production activity 

may involve several reactors connected with inflow and outflow pipes. If there are n reactors, 

the whole process is governed by a system of 𝑛 differential equations of the form 

[

𝑥′(𝑡)

𝑦′(𝑡)

𝑧′(𝑡)
] =

[
 
 
 
 −

8

3
−
4

3
1

−
17

3
−
4

3
1

−
35

3

14

3
−2]
 
 
 
 

[

𝑥1
𝑥2
𝑥3
] + [

12
29
48
]  

Analytical Solution 



Int. J. Anal. Appl. 19 (6) (2021) 942 

 

𝑥(𝑡) =
1

6
𝑒−3𝑡(6 − 50𝑒𝑡 + 10𝑒2𝑡 + 34𝑒3𝑡),  

𝑦(𝑡) =
1

6
𝑒−3𝑡(12 − 125𝑒𝑡 + 40𝑒2𝑡 + 73𝑒3𝑡),  

𝑧(𝑡) =
1

6
𝑒−3𝑡(14 − 200𝑒𝑡 + 70𝑒2𝑡 + 116𝑒3𝑡). See [5] for more info. 

 

Application Problem 2 

We examine the coefficient matrix A with the initial conditions 𝑦(0) and the forcing function 

𝑔(𝑡) given by 

𝐴 = [
0 −1 1
0 2 0
−2 −1 3

], 

𝑦(0) = [
0
0
1
], 

𝑔(𝑡) = [
1
𝑡
𝑒−𝑡

]. 

Analytical Solution 

𝑦(𝑡) =

[
 
 
 
 
7

12
𝑒2𝑡 +

3

2
𝑒𝑡 +

𝑒−𝑡

6
−

𝑡

2
−
9

4

𝑒2𝑡

4
−

𝑡

2
−
1

4

17

12
𝑒2𝑡 +

3

2
−
𝑒−𝑡

6
−

𝑡

2
−
7

4 ]
 
 
 
 

. 

 See [28] for more info. 

 

Application Problem 3 

We study the initial value problem 𝑦′(𝑡) = 𝐴𝑦(𝑡),𝑦(0) = [1,0,−1]𝑇, where 

 

𝐴 = [
−21 19 −20
19 −21 20
40 −40 −40

]. 

Analytical Solution 

𝑦(𝑡) =

[
 
 
 
1

2
𝑒−2𝑡 +

1

2
𝑒−40𝑡(𝑐𝑜𝑠(40𝑡) + sin(40𝑡)))

1

2
𝑒−2𝑡 −

1

2
𝑒−40𝑡(𝑐𝑜𝑠(40𝑡) + sin(40𝑡)))

−𝑒−40𝑡 + (𝑐𝑜𝑠(40𝑡) − sin(40𝑡))) ]
 
 
 
.  

See [18] for more info. 

 



Int. J. Anal. Appl. 19 (6) (2021) 943 

 

IV.  RESULTS 

This aspect implements the computational strategy of variable step, variable order for 

solving stiff systems of ordinary differential equations. The different application problems of 

stiff system of ordinary differential equations were implemented on diverse convergence 

criteria of the following: 10−3,10−3,10−4,10−5,10−6,10−7,10−8,10−9,10−10 𝑎𝑛𝑑 10−11 . The 

results of the application problems 1, 2 and 3 are displayed on tables 1, 2 and 3.  

 

TABLE I. CSVSVO IMPLEMENTATION 

 

𝑀𝑒𝑚𝑝𝑙𝑜𝑦𝑒𝑑 𝑀𝑎𝑥𝑒𝑟𝑟𝑜𝑟𝑠 𝐶𝑜𝑛𝑣𝑐𝑟𝑖𝑡𝑒𝑟𝑖𝑎 

CSVSVO 6.73296E-05 10−2 

CSVSVO 1.50195E-01 10−2 

CSVSVO 2.30889E-01 10−2 

CSVSVO 7.17852E-09 10−4 

CSVSVO 1.81011E-03 10−4 

CSVSVO 3.07211E-03 10−4 

CSVSVO 7.17852E-13 10−6 

CSVSVO 1.84418E-05 10−6 

CSVSVO 3.16097E-05 10−6 

CSVSVO 6.24284E-16 10−8 

CSVSVO 1.84762E-07 10−8 

CSVSVO 3.17E-07 10−8 

CSVSVO 1.96457E-16 10−10 

CSVSVO 1.84796E-09 10−10 

CSVSVO 3.17E-09 10−10 

 

 

 

 

 

 



Int. J. Anal. Appl. 19 (6) (2021) 944 

 

 

TABLE II. CSVSVO IMPLEMENTATION 

 

𝑀𝑒𝑚𝑝𝑙𝑜𝑦𝑒𝑑 𝑀𝑎𝑥𝑒𝑟𝑟𝑜𝑟𝑠 𝐶𝑜𝑛𝑣𝑐𝑟𝑖𝑡𝑒𝑟𝑖𝑎 

CSVSVO 1.33456E-04 10−3 

CSVSVO 1.10578E-02 10−3 

CSVSVO 3.81005E-02 10−3 

CSVSVO 1.43869E-08 10−5 

CSVSVO 2.10822E-04 10−5 

CSVSVO 4.47082E-04 10−5 

CSVSVO 1.43996E-12 10−7 

CSVSVO 2.22913E-06 10−7 

CSVSVO 4.47991E-06 10−7 

CSVSVO 1.11022E-16 10−9 

CSVSVO 2.24143E-08 10−9 

CSVSVO 4.48E-08 10−9 

CSVSVO 1.11022E-16 10−11 

CSVSVO 2.24226E-10 10−11 

CSVSVO 4.48E-10 10−11 

 

 

 

 

 

 

 

 

 

 

 

 



Int. J. Anal. Appl. 19 (6) (2021) 945 

 

 

TABLE III. CSVSVO IMPLEMENTATION 

 

𝑀𝑒𝑚𝑝𝑙𝑜𝑦𝑒𝑑 𝑀𝑎𝑥𝑒𝑟𝑟𝑜𝑟𝑠 𝐶𝑜𝑛𝑣𝑐𝑟𝑖𝑡𝑒𝑟𝑖𝑎 

CSVSVO 6.73296E-05 10−2 

CSVSVO 1.50195E-01 10−2 

CSVSVO 2.30889E-01 10−2 

CSVSVO 7.17852E-09 10−4 

CSVSVO 1.81011E-03 10−4 

CSVSVO 3.07211E-03 10−4 

CSVSVO 7.17852E-13 10−6 

CSVSVO 1.84418E-05 10−6 

CSVSVO 3.16097E-05 10−6 

CSVSVO 6.24284E-16 10−8 

CSVSVO 1.84762E-07 10−8 

CSVSVO 3.17E-07 10−8 

CSVSVO 1.96457E-16 10−10 

CSVSVO 1.84796E-09 10−10 

CSVSVO 3.17E-09 10−10 

 

 

V. CONCLUSION 

Applications problem 1, 2 and 3 represents Stiff systems of ordinary differential 

equations which seems to generates an unstable system behavior, and as such requires a 

technical approach like CSVSVO to guarantee an improve efficiency and better accuracy. The 

stiff systems of ordinary differential equations are carried out employing the CSVSVO 

implementation. The CSVSVO has the capacity to introduce the convergence criteria in order to 

engender the desired result is achieved. These convergence criteria decide whether the result is 

accepted or rejected. The CSVSVO implementation is done utilizing a multiprocessor approach 

executed under the Mathematica Software platform. Tables 1, 2 and 3 displays the 



Int. J. Anal. Appl. 19 (6) (2021) 946 

 

computational results established that the CSVSVO is reached via the convergence criteria of 

the following: 10−3,10−3,10−4,10−5,10−6,10−7,10−8,10−9 

,10−10 𝑎𝑛𝑑  10−11. In addition, with the trend of the maximum errors achieved via the different 

convergence criteria, we can conclude that the CSVSVO is capable of resolving stiff systems of 

ordinary differential equations with better efficiency and accuracy as exhibited in Tables 1, 2 

and 3.  See [5, 18, 21—25, 28] for details.  

 

Acknowledgements: The authors would like to appreciate Covenant University for providing 

sponsorship throughout the study period of time. Many thanks to the anonymous reviewers for 

their continuous contribution. 

 

Conflicts of Interest: The author(s) declare that there are no conflicts of interest regarding the 

publication of this paper. 

 

REFERENCES 

[1] M. L. Abell, J. P. Braselton, Differential Equations with Mathematica, Elsevier Academic 

Press, USA, 2004 

[2] U. M. Ascher, L. R. Petzoid, Computer Methods for Ordinary Differential Equations and 

Differential-Algebraic Equations, SIAM, USA, 1998.  

[3] K. Atkinson, W. Han, D. Stewart, Numerical solution of ordinary differential equations, 

John Wiley & Sons, Inc., New Jersey, 2009. 

[4] F. Ceschino, Modification de la longueur du pas dans l’ integration numerique par les 

methods a pas lies, Chifres, Vol. 2, 101-106, 1961. 

[5] W. Cheney, D. Kingaid, Numerical Mathematics and Computing, Thomson Brooks/Cole, 

USA, 2008. 

[6] J. R. Dormand, Numerical Methods for Differential Equations, CRC Press, New York, 1996.  

[7] J. D. Faires, R. L. Burden, Initial-Value Problems for ODEs, Dublin City University, 2012.  

[8] S. O. Fatunla, Numerical Methods for Initial Value Problems in Ordinary Differential 

Equations, Academic Press, Inc., New York, 1988. 

[9] C. V. D. Forrington, Extensions of the predictor-corrector method for the solution of 

systems of ordinary differential equations, Comput. J. 4 (1961), 80-84. 



Int. J. Anal. Appl. 19 (6) (2021) 947 

 

[10] C. W. Gear, Numerical Value Problems in ODEs, Prentice-Hall, Inc., New Jersey, USA, 

1971. 

[11] C. W. Gear, K. W. Tu, The effect of variable mesh size on the stability of multistep methods, 

SIAM J. Numer. Anal. 11 (1974), 1025-1043. 

[12] C. W. Gear, D. S. Watanabe, Stability and convergence of variable order multistep methods, 

SIAM J. Numer. Anal. 11 (1974), 1044-1058. 

[13] E. Hairer, S. P. Norsett, G. Wanner, Solving ordinary differential equations I, Springer 

Heidelberg Dordrecht, New York, 2009. 

[14] M. I. Hazizah, B. I. Zarina, Diagonally implicit block backward differentiation formula with 

optimal stability properties for stiff ordinary differential equations, Symmetry, 11 (2019), 

1342. 

[15] F. T. Krogh, A variable step variable order multistep method for the numerical solution of 

ordinary differential equations, Information Processing 68 North-Holland, Amsterdam, 

194-199, 1969. 

[16] F. T. Krogh, Algorithm for changing the step size, SIAM J. Num. Anal. 10 (1973), 949-965. 

[17] F. T. Krogh, Changing step size inthe integration of differential equations using modified 

divided differences, Proceedings of the Conference on the Num. Sol. of ODE, Lecture Notes 

in Math Vol. 362, 22-71, 1974. 

[18] J. D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley & 

Sons, Inc., New York, 1973. 

[19] J. D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, 

Inc., New York, 1991. 

[20] C. B. Moler, Numerical Computing with MATLAB, SIAM, USA, 2004. 

[21] J. G. Oghonyon, S. A. Okunuga, N. A. Omoregbe, O. O. Agboola, A computational 

approach in estimating the amount of pond and determining the long time behavioural 

representation of pond pollution”, Global Journal of Pure and Applied Mathematics, Vol. 

11, 2773-2786, 2015. 

[22] J. G. Oghonyon, J. Ehigie, S. K. Eke, Investigating the convergence of some selected 

properties on block predictor-corrector methods and its applications, J. Eng. Appl. Sci. 11 

(2017), 2402-2408. 



Int. J. Anal. Appl. 19 (6) (2021) 948 

 

[23] J. G. Oghonyon, O. A. Adesanya, H. Akewe, H. I. Okagbue, Softcode of multi-processing 

Milne’s device for estimating first-order ordinary differential equations, Asian J. Sci. Res. 11 

(2018), 553-559. 

[24] J. G. Oghonyon, O. F. Imaga, P. O. Ogunniyi, The reversed estimation of variable step size 

implementation for solving nonstiff ordinary differential equations, Int. J. Civil Eng. 

Technol. 9 (2018), 332-340. 

[25] J. G. Oghonyon, S. A. Okunuga, H. I. Okagbue, Expanded trigonometrically matched block 

variable-step-size technics for computing oscillating vibrations, Lecture Notes in 

Engineering and Computer Science, Vol. 2239, 552-557, 2019. 

[26] P. Piotrowsky, Stability, consistency and convergence of variable k-step methods for 

numerical integration of large systems of ordinary differential equations, Lecture Notes in 

Math., Vol. 109, 221—227, 1969. 

[27] S. E. Ramin, Numerical Methods for Engineers and Scientists Using MATLAB®, CRC Press, 

London, 2017. 

[28] P. M. Shankar, Differential Equations: A problem solving approach based on MATLAB®, 

CRC Press, USA, 2018.