Numerical solutions for ordinary differential equations (ODEs) are very important in scientific computation, as they are widely used to model real world problems      . Indeed, ODEs are found in studies of electrical circuits, chemical kinetics    , vibrations, simple pendulum  , rigid body rotation  , atmospheric chemistry problems    , motions of the planet in a gravity field like the Kepler problem   , biosciences and many more fields. However, most realistic models developed from these equations cannot be solved analytically, especially in case of nonlinear differential equations. The study of such models therefore requires the use of numerical methods. Nevertheless, there are several ODEs, known as “stiff” equations, which classical numerical methods do not handle very efficiently.
The phenomenon of stiffness was first recognized by Curtiss and Hirschfelder  in 1952. Since then, an enormous amount of effort has gone into the analysis of stiff problems and, as a result, a great many numerical methods have been proposed for their solution.
Stiff problems are too important to ignore, and are too expensive to overpower. They are too important to ignore because they occur in many physically important situations. They are too expensive to overpower because of their size and the inherent difficulty they present to classical methods. Even if one can bear the expense, classical methods of solution require so many steps that round off errors may invalidate the solution  .
Stiff problem entails rapidly decaying transient solution, which arises naturally in wide variety of applications including the study of spring and damping systems, the analysis of control system and problems in the chemical kinetics   . Stiff differential equations also occur in other kind of studies, such as biochemistry, biomedical systems, weather prediction, mathematical biology and electronics. The importance of stiff equations is discussed by Shampine and Gear  and Bjurel et al.  , who present a comprehensive survey of application areas in which stiff equations arise.
We have to emphasize that while the intuitive meaning of the term stiff is clear to all specialists, much controversy is going on about its mathematical definition. The main purpose of this paper is to outline some of the important characteristics of stiff problems and to present the explicit one-step algorithm suggested by Fatunla  for numerical integration of stiff and highly oscillatory differential equations, with an application to the Robertson problem (RP)   . To the best of our knowledge there is no mention of this application in the literature. Note that the RP consists of a system of three non-linear ODEs modeling the kinetics of three species. It is very popular in numerical studies     and is often used as a test problem in the stiff integration comparisons.
The Explicit Fatunla’s method (EFM) has been used by Frapiccini et al.  to solve numerically the time-dependent Schrödinger equation describing physical processes whose complexity requires the use of state-of-the-art methods.
The rest of this paper is organized as follows. In Section 2, we concentrate on some basic properties of stiff differential equations. Section 3 contains a brief introduction to the EFM, which has been shown to be available for solving stiff ODEs   . We show how to use this method when one is confronted with first-order stiff systems of scalar ODEs. In Section 4, we apply the method in question to the RP and compare the obtained results with those given by the solver RADAU  which is based on implicit Runge-Kutta methods. Finally, we present our conclusions in Section 5.
2. Stiff Differential Equations and Stability Analysis
As stated above, stiff differential equations appeared in the early 1950s as a result of some pioneering work by Curtiss and Hirschfelder  , and some ten years later, one could say, in the words of Dahlquist   , that “… Around 1960, things became completely different and everyone became aware that the world was full of stiff problems”.
Stiffness is a subtle, difficult, and important concept in the numerical solution of ODEs. It depends on the differential equation, the initial conditions, and the numerical methods  . Dictionary definitions of the word stiff involve terms like “not easily bent”    , “rigid”   , “stubborn”    and “hard”  .
The differential equations we consider in this work are of the form
where x is the independent variable which often plays the role of time (t), and is an unknown function that is being sought. The given function of two variables defines the differential equation. This equation is called a first-order differential equation because it contains a first-order derivative. Sometimes, for convenience, we will omit the x argument in .
Let us indicate that the general solution of the first-order Equation (1) normally depends on an arbitrary integration constant. To single out a particular solution, we need to specify an additional condition. Usually such a condition is taken to be of the form
The differential Equation (1) and the initial value condition (2) together form an initial value problem (IVP):
We have to emphasize that many physical applications lead to higher-order systems of ODEs, but there is a simple reformulation that can convert them into equivalent first order systems  . Thus, we do not lose any generality by restricting our attention to the first-order case throughout this section.
2.1. Stability Analysis
Examining the stability question for the general problem (3) is too complicated. Instead, we examine the stability of numerical methods for the model problem
whose analytical solution is
where is a constant. The equation
is called the “Dahlquist test equation”  . The simple test problem (4) has traditionally been used for stability analysis since we can easily obtain analytical expressions describing the solution produced by the numerical method. Studying the behavior of a numerical method in solving this problem is also useful in predicting its behavior in solving the general problem (3). Indeed, if we expand about , we obtain
with and .
This is a valid approximation if is sufficiently small, i.e., where h is a small positive real number. Introducing , we obtain
The inhomogeneous term will drop out of all derivations concerning numerical stability, because we are concerned with differences of solutions of the equation. Dropping in Equation (9), we obtain the model Equation (6).
If we are dealing with a set of m equations, i.e.
where is a vector of m components, and is a vector-valued function of variables, the partial derivative intervening in Equation (7) becomes a matrix called the Jacobian of . Thus the model equation becomes
where is the Jacobian of and can be regarded as a constant matrix. For a linear system of differential equations, , determined by the matrix , the Jacobian is precisely the matrix . The differential system (11) reduces to a set of m equations like (6), i.e.
with the eigenvalues of and , , the components of the vector defined by
where is a matrix of right eigenvectors of the Jacobian, which means that .
We now study the behavior of the solution computed by the explicit Euler method (EEM)   for the test problem (4) using a fixed step size h. We assume that we have computed solution values at points , respectively, where for . The above mentioned numerical method, when applied to the model problem (4), leads to the following equation:
The ratio of the computed solutions at and is given by
We compare this with the ratio of the true solutions at the same points, which is
It is clear that when is real, is a reasonable approximation to except if . For large negative , is much smaller than 1, while is (in magnitude) greater than 1. This means that the numerical solution is growing while the true solution is decaying. We say that the numerical solution produced by the EEM for is unstable.
Another way of looking at stability is to look at the propagation of local errors. For the EEM we have    :
where . If , we say that the errors are growing if . In other words, the errors grow if and, for such step sizes, the method is called unstable.
To sum up, the requirement as , which, for , mirrors the asymptotic behavior of the analytical solution (5), implies that the inequality
must be satisfied. If , then the step size of integration is severely restricted by Equation (18) even though the true solution becomes negligible very quickly.
Before finishing this subsection, let us define some stability concepts for one-step numerical methods. For this purpose, we consider the model problem (4) where is a constant (possibly complex) and . If we assume that , the exact solution of this problem is
It is clear that if and only if . We have therefore to look for the conditions that have to be imposed on the above mentioned numerical methods in order that the numerical solution as where h is the step size of integration. By applying any one-step numerical method to the model problem in question, we obtain  
where and is the so-called stability function. In order that tends to zero as , we must impose thereby implying some constraints on the step size h. The set is called the stability domain of the numerical method. This latter one is said to be A-stable (absolutely stable)   if its stability domain is included in . It is L-stable  if, apart from being A-stable, the stability function has the property . A numerical scheme is called -stable   for if its region of absolute stability includes the infinite wedge
This definition calls for a numerical method to have a stability region which is unbounded but which does not include the whole complex left hand half-plane. One of the reasons why this is such a useful definition is that many problems have eigenvalues of the Jacobian that lie in a sector .
L-stable methods are the most stable ones  and the backward Euler method   is an example of an A-stable method.
We now briefly discuss stiffness and the difficulties involved in solving stiff equations. As discussed earlier, a system of ODEs of the form may be approximated by over a small interval of the independent variable x. If is diagonalizable, this new system may be transformed into
where is a diagonal matrix with eigenvalues of (possibly complex) on its diagonal, and is related to by Equation (13).
Many differential equation systems of practical importance in scientific modeling exhibit a distressing behavior when solved by classical numerical methods. This behavior is distressing because these systems, classified as stiff, are characterized by very high instability when approximated by standard numerical methods.
An intuitive idea of what stiff equations are is that they are problems with some smooth and some transient solutions, where all solutions reach the smooth one after a short time (after the transient phase has finished)  .
Since stiffness is closely related to the behavior of perturbations to a given solution, it is important to study the effect of a small perturbation to a solution . The parameter is small, in the sense that we are interested only in asymptotic behavior of the perturbed solution as this quantity approaches zero. If is replaced by in the differential equation and the solution expanded in a series in power of , with and higher powers replaced by zero, we obtain the system
where is the Jacobian of , and . Subtracting Equation (10) from Equation (23), we obtain the equation
which controls the behavior of the perturbation. The Jacobian matrix has a crucial role in the understanding of stiff problems. In fact, its spectrum is sometimes used to characterize stiffness    . In an interval , chosen so that there is a moderate change in the value of the solution to , and very little change in , the eigenvalues of determine the growth rate of components of the perturbation. The existence of one or more large and negative values of , for , the spectrum of , indicates that stiffness is almost certainly present. If possesses complex eigenvalues, then we interpret this test for stiffness as the existence of a such that is negative with large magnitude. A definition of stiffness in terms of the spectrum of the Jacobian matrix has been proposed by Gaffney   who, in 1984 discussed the concept of stiffness oscillatory systems. This definition states that the IVP
is considered to be stiff oscillatory over the interval S if for every the eigenvalues of the Jacobian possess the following properties:
or if the stiffness ratio satisfies
and for at least one pair of j in .
If an explicit method like the Euler’s one is used to solve the differential system when has a real eigenvalue , the step size of integration is controlled by a transient solution (which quickly becomes negligible), whereas outside the transient phase we wish the step size to be controlled by accuracy alone. This suggests a rather more pragmatic definition of stiffness, where the definition is not based on an analysis of the actual differential equation to be solved but is instead based on the relative performance of implicit and explicit integration methods. Perhaps the best and the oldest definition of this type is due to Curtiss and Hirschfelder  who, in 1952, said: “Stiff equations are equations where certain implicit methods, and in particular backward-differentiation formulae (BDFs), perform better, usually tremendously better, than explicit methods”.
The same situation can occur for complex eigenvalues with negative real parts. The problem is that the EEM is not A-stable or even -stable for any . The same is true for all explicit methods like Euler’s one: no such explicit method can be -stable. We are therefore forced to use implicit methods, like the backward Euler method, to solve stiff systems. These methods require more work per step than explicit methods, however, since a system of nonlinear algebraic equations must be solved at each step. For further details of stiff problems, the reader is referred to Shampine and Gear  , Lambert  and Söderling et al.  .
Before finishing this section, let us introduce a definition of stiffness which involves a certain norm that depends on the differential system to be solved. If , the Jacobian matrix of the system, is diagonalizable, this definition is as follows. A system of ODEs is stiff if its stiff indicator, defined by
is “large” and negative. Here, and are the greatest lower bounds (glb) and the least upper bound (lub) logarithmic Lipschitz constants, respectively   . Note that corresponds to the usual logarithmic norm  .
The stiffness indicator is easily computed. Let denote the eigenvalues of a matrix . Then, for the Euclidian norm, it holds that
where denotes the hermitian part of the matrix . Other norms can also be used, provided that and are computed using well-known expressions.
3. Explicit Fatunla’s Method
The idea behind the EFM is to take into account the stiffness of the differential system by introducing the stiffness parameters in interpolating functions that approximate the solution. This allows one to deal with differential systems where the Jacobian matrix displays eigenvalues (possibly complex with large negative real part) that differ by many orders of magnitude. That explains why Fatunla’s method has the capacity to solve stiff equations.
In his paper  , Fatunla considers initial value problems of type
with in the finite interval , where for some positive integer . It is assumed that is sufficiently differentiable. Throughout this section, we use the same vector notation as in reference  . More precisely, we write and . The numerical estimates to the theoretical solution at the points , are to be generated.
According to the EFM, the theoretical solution is, on every subinterval , approximated by the interpolating function
and being m-tuples with real entries, is the identity matrix, whilst and are diagonal matrices, usually called the stiffness matrices; or
where and are m-tuples with complex entries, and (*) denotes complex conjugate. The choice of interpolation is determined by Equation (39).
3.1. Case of Real Interpolation Formula
By demanding that the interpolating function (32) coincides with the theoretical solution at the endpoints of the interval , the following recursion formula is readily obtainable    :
where we use the notations , . and represent diagonal matrices defined by
where and are diagonal matrices with non zero entries in the main diagonal given by
The components of the stiffness matrices are given by   :
where and , ( ) are expressed in terms of the components of the derivatives , :
provided that the denominator in Equation (38) is not zero.
3.2. Case of the Complex Interpolation Formula
The complex interpolation formula (see Equation (33)) is adopted if in Equation (37), the following relationship holds:
The components of the stiffness matrices, which, in this case, are and , can be written as
for some real numbers and . By imposing the same restrictions on (33) as (32), we still obtain the interpolation formula (34) but now with components and of and , respectively, given by 
3.3. Local Truncation Error
Fatunla  has shown that his method, when applied to the scalar test problem , is L-stable and exponentially fitted for any complex value with negative real part. This means that the corresponding stability function gives the optimum stability properties. Furthermore, it can be shown that the jth component of the local truncation error at is given by   :
for the real interpolation formula and by 
for the complex interpolation formula.
4. Application to the Robertson Problem
As stated earlier, the RP consists of a system of three non-linear ODEs, i.e.
modeling the kinetics of three species, with and parameters , and . , and denote the concentrations of the three species. The RP is well-known to be stiff    . Note that means . The Jacobian matrix is here given by
Its Hermitian part is
This one has certainly three real eigenvalues which are solution of the cubic equation
where , , and
Using the Cardan method for cubic equations in one variable, we obtain
, , (49)
with , and , where p and q are defined by , . We can therefore study the stiffness of the RP by use of the stiffness indicator (see Equation (29))
where and are, respectively, the largest and smallest eigenvalues of .
The components of the solution obtained by means of the EFM is shown in Figure 1. This figure contains also plots of the same components computed by the RADAU solver  which is based on implicit Runge-Kutta methods. In Figure 2, we show the stiffness indicator , together with the variation of the integration step size associated with the EFM.
We want to emphasize that implementation of the EFM to solve the RP requires the calculation of the function and its four first derivatives, i.e. , , and at each value . The function is here given by
We have used maple 18 software to establish analytical expression of , , and .
As expected, we remark, from the plot of the stiffness indicator, that the RP is very stiff. We also find that our results agree noticeably with those computed by the RADAU solver except for the solution component for small values of the time. We think that the observed discrepancy between our results and the
Figure 1. The Robertson equation. Solution components y1 and y3 (top) and the component y2 (bottom).
Figure 2. Stiffness indicator (top) and the integration step size (bottom) for the EFM.
RADAU ones for t at around 10−4 atomic units (a.u.) is due to the fact that during the implementation of the EFM, we have chosen, for a.u., a very small truncation error (10−12 in magnitude) to control the integration step size, which means that our results are probably more accurate. In the case of a.u., the time step has been adapted according to the condition . It is because we have used two distinct conditions on the truncation error for the control of the integration step size that we have a discontinuity at a.u. in the graph of the step size.
The object of this paper has been to present some basic characteristics of stiff differential equations, as well as to introduce the EFM which has been shown to be available for solving stiff problems. We have shown that the stiffness indicator of the Jacobian matrix gives a sufficient information to estimate the computational costs of explicit schemes like the Euler’s one. Numerical results obtained solving the Robertson problem using the EFM on the one hand and the solver RADAU on the other hand, confirm the fact that the explicit method in question can be a good candidate to solve stiff ODEs.
We have to emphasize that modern codes for solving ODEs automatically vary the step size, estimate the local error, and provide facilities to compute the solution at intermediate points via interpolation  . That is why during the implementation of the EFM we calculate the truncation error to control the size of the integration step imposing a boundary criterion for .
 Mazzia, F. and Iavernaro, F. (2003) Test Set for Initial Value Problem Solvers, Release 2.4. http://www.dm.uniba.it/~testset
 Akinfenwa, O., Jator, S. and Yoa, N. (2011) An Eight Order Backward Differential Formula with Continuous Coefficients of Stiff Ordinary Differential Equations. International Journal of Mathematical and Computer Sciences, 7, 171-176.
 Curtiss, C.F. and Hirschfelder, J.O. (1952) Integration of Stiff Equations. Proceedings of the National Academy of Sciences of the United States of America, 38, 235-243. https://doi.org/10.1073/pnas.38.3.235
 Sandu, A., Verwer, J.G., Blom, J.G., Spee, E.J., Carmichael, G.R. and Potra, F.A. (1997) Benchmarking Stiff ODE Solvers for Atmospheric Chemistry Problems II: Rosenbrock Solvers. Atmospheric Environment, 31, 3459-3479.
 Kin, J. and Cho, S.Y. (1997) Computational Accuracy and Efficiency of the Time-Splitting Method in Solving Atmospheric Transport/Chemistry Equations. Atmospheric Environment, 31, 2215-2224.
 Ebady, A.M.N., Habib, H.M. and El-Zahar, E.R. (2012) A Fourth Order A-Stable Explicit One-Step Method for Solving Stiff Differential Systems Arising in Chemical Reactions. International Journal of Pure and Applied Mathematics, 81, 803-812.
 Hairer, E. and Wanner, G. (1999) Stiff Differential Equations Solved by Randau Methods. Journal of Computational and Applied Mathematics, 111, 93-111.
 Anake, T.A., Bishop, S.A., Adesanya, A.O. and Agarana, M.C. (2014) An -Stable Method for Solving Initial-Value Problems of Ordinary Differential Equations. Advances in Differential Equations and Control Processes, 13, 21-35.
 Cash, J.R. (2003) Review Paper: Efficient Numerical Methods for the Solution of Stiff Initial-Value problems and Differential Algebraic Equations. Proceedings of the Royal Society A, 459, 797-815. https://doi.org/10.1098/rspa.2003.1130
 Gaffney, P.W. (1984) A Performance Evaluation of Some FORTRAN Subroutines for the Solution of Stiff Oscillatory ODEs. ACM Transactions on Mathematical Softwares, 10, 58-72. https://doi.org/10.1145/356068.356073