INTRODUCTION
Generally speaking, a special third order differential equations (ODEs) of the form:
with initial conditions:
y(t_{0}) = α, y'(t_{0}) = β, y^{''}(t_{0}) = γ
where, f:RxR→R which is not explicitly dependent on the first derivative y'(x) and the second derivative y''(x) of the solution. The ordinary differential Eq. 1 is frequently found in many physical problems such as thin film flow, gravity and electromagnetic waves. Most researchers, scientists and engineers used to solve Eq. 1 by converting the third order differential equations to a system of first order equations three times the dimension (Mechee et al., 2013a). However, it is more efficient if the problem can be directly solved using numerical methods. Such a type of work can be seen in Awoyemi and Idowu (2005), Waeleh et al. (2011), Zainuddin (2011) and Jator (2010). All methods previously discussed are multistep methods; hence they need the starting values when used to solve ODEs Eq. 1. Most of the methods for solving ODEs can also be adapted for solving delay differential equations (DDEs).
In recent years there has been a growing interest in numerical solutions of DDEs, this is due to the appearance of such equations in various areas such as neural network theory, epidemiology and time lag control processes (Mechee et al., 2013b). DDEs also provide us with realistic model of many phenomena arising in real world problems for example DDEs can be used in modeling of population dynamics and spread of infectious diseases and two body problems of electrodynamics (Bellen and Zennaro, 2003; Forde, 2005; Driver, 1977; Smith, 2011; Erneux, 2009; Kuang, 1993). In this study, we are concerned with the onestep method particularly the RungeKutta method of order three for directly solving third order ordinary differential equations. Accordingly, we have developed a direct RungeKutta (RKD) method which can be directly used to solve Eq. 1. The advantage of the new method over multistep methods is that it is self starting. The method is then adapted for solving special third order DDEs with multiple delays. The third order DDE can be written in the following form:
with initial conditions:
where, φ is a continuous function and τ_{1}, τ_{2}, …, τ_{n }are time delays.
Numerical results on two sets of problems consisting of ordinary and delay differential equations are given and compared with the numerical results when the problems are reduced to a system of first order ODEs and DDEs, respectively and solve using RungeKutta methods. Stability polynomial of the method when applied to linear third order DDE is also presented. Numerical results on two sets of problems consisting of ordinary and delay differential equations are given and compared with the numerical results when the problems are reduced to a system of first order ODEs and DDEs respectively and solve using RungeKutta methods.
Derivation of RKD method
The general form of RKD method with sstage for solving initial value problem Eq. 1 can be written as:
Where:
for i = 2, 3, ..., s. The parameters of RKD method are for i = 1, 2, ..., s and j = 1, 2, ..., s are assumed to be real. If a_{ij} = 0 for i<j, it is an explicit method and otherwise implicit method. The RKD method can be expressed in Butcher notation using the table of coefficients as follows:
To determine the coefficients of the RKD method, the expressions given in Eq. 37 are expanded using Taylor's series expansion. After some algebraic manipulations this expansion a equated to the true solution which are given by Taylor's series expansion. General order conditions for the RKD method can be found from the direct expansion of the local truncation error. The order conditions can found in Mechee et al. (2013a) which introduced direct RungeKutta (RKD) method of order five with three stage for solving thin film problem.
ORDER CONDITIONS OF THE METHOD
Mechee et al. (2013b) derived the order condition of RKD method up to order six. In this study using the same technique to get the order conditions up to order four. The order conditions for twostage third order RKD method can be written as follows.
Order conditions for y:
Order conditions for y':
Order conditions for y'':
All indices are from 1 to s. To get thirdorder RKD method, the following simplifying assumption is used in order to reduce the number of equations to be solved:
ZERO STABILITY OF THE METHODS
Next, we will discuss the zerostability of the methods it is one of the criteria for the method to be convergent. Zerostability is an important tool for proving the stability and convergence of linear multistep methods. The interested reader is referred to the textbooks by Lambert (1991) and Butcher (2008) in which zerostability is discussed. Zerostability has been discussed in Hairer et al. (2010), where it is used to determine an upper bound on the order of convergence of linear multistep methods.
In studying the zero stability of RKD method, we can write method, Eq. 24 as follows:
Thus the characteristic polynomial is:
Hence, the method is zerostable since the roots are ∈ = 1, 1, 1, are less or equal to one.
DERIVATION OF THIRD ORDER RKD METHODS
The RKD method of sstage and porder can be derived by the solving the order conditions of the method. The system of nonlinear equations (order conditions) of the method depend on p, however, the existence of the solutions of this system depends on the number of coefficients of the method which depends on the sstage of the method in addition to the number of independent order conditions of the method.
DERIVATION OF TWOSTAGE THIRDORDER RKD METHODS
To derive the twostage and thirdorder RKD method, we use the algebraic conditions up to order three in the equations of order conditions in y, y' and y'' in Eq. 8, 10, 11 and 1315. The resulting system of equations consists of six nonlinear equations with 8 unknowns variables to be solved. Consequently, there is a solution with two free parameters b_{1}:
and a_{21} however we chose one of them arbitrary;
but the second free coefficients can be chosen using minimization of the truncation error. Accordingly Dormand (1996) the free parameters can be chosen by minimizing the global error of the fourth order conditions minimize the global error. The technique is as follows:
• 
First: We find the error coefficients of y, y' and y'' respectively as the following: 
Table 1:  Butcher tableau 

• 
Second: We find the global error norm as the following: 
Finally, we minimize the four truncation errors in Eq. 2124 with respect to the free parameter b_{1}, in this case we get the free parameter as b_{1} = 1/10. The error norms for y_{n}, y'_{n} and y''_{n} are givern by:
respectively, where τ^{(4)}, τ’^{(4)} and τ’=^{(4)} are error terms of the fourthorder conditions for y, y' and y'' respectively. The RKD method of order three and twostage is denoted by RKD3 which can be expressed in the following Butcher tableau in the Table 1.
Adapting RKD methods for directly solving third order DDEs: Consider problem Eq. 2 for initial value problem for delay differential equations. To find the solution y_{i+1} at the point t_{i+1} for i = 0, 1, …, n1, we need the values of solutions at point t_{i} and all time delays points t_{i}τ_{1}, t_{i}τ_{2}, t_{i}τ_{3}, …, t_{i}τ_{n} for i = 0,1,...,n (Mechee et al., 2013a). The general form of RKD methods with sstage for solving the initial value problem DDEs Eq. 2 can be written as the following:
The parameters of the methods Eq. 2529 are and a_{ij } for i = 1, 2, 3, …, s and j = 1, 2, 3, …, s are given in Table 1.
TIME DELAY INTERPOLATION
Let the interval of the differential Eq. 2 be I = [a, b] and t_{i} = a+ih for i = 0, 1, 2, …, n and:
where n the number of points in the interval I.
The numerical method approximates the solution y_{i+1} at the point t_{i+1} for i = 0, 1, 2, …, n1,
Hence:
To evaluate the delay term we used cubic interpolation and the details of it can be referred to Mechee et al. (2013b).
STABILITY OF THE METHODS WHEN APPLIED TO THIRD ORDER DDE
To study the stability of numerical method (Eq. 35) , consider the linear test equation:
when, the method is applied to the linear test Eq. 31, (Jiaoxun and Yuhao, 2005; Hongjiong and Jiaoxun, 1995; AlMutib, 1984; Liu and Spijker, 1990), we have:
Where:
and:
such that:
So:
and:
where, H_{α} = (αh)^{3}, H_{μ} = (μh)^{3}.
So:
where:
Hence, the stability polynomial of the method is:
where:
NUMERICAL RESULTS
In this section a set of third order ordinary and delay differential equations are solved using RKD3 method of order three and numerical results are compared with the existing RK methods of the orders three and four.
The following notations are used in Fig. 18:
• 
h: Stepsize used 
• 
RKD3: The Direct RungKutta method of order three derived in section 5 
• 
RK3: Existing RungeKutta method order three as given in Butcher (2008) 
• 
RK4: Existing RungeKutta method order fourth as in Dormand (1996) 
• 
Total time: The total time in second to solve the problems 
• 
MAX ERROR: Max y(x_{n})y_{n} Absolute value of the true solution minus the computed solution 

Fig. 1:  A log max of errors versus computational time for problem 1 

Fig. 2:  A log max of errors versus computational time for problem 2 

Fig. 3:  A log max of errors versus computational time for problem 3 

Fig. 4:  A log max of errors versus computational time for problem 4 

Fig. 5:  A log max of errors versus computational time for problem 5 

Fig. 6:  A log max of errors versus computational time for problem 6 

Fig. 7:  A log max of errors versus computational time for problem 7 

Fig. 8:  A log max of errors versus computational time for problem 8 
PROBLEMS OF ODES
Problem 1(Homogenous linear):
y'''(t) = y(t), 0<t<b
Y(0) = 1, y'(0) = 1, y' '(0) = 1,
y(t) = e^{t}, b = 1
Problem 2(Non homogenous linear):
y'''(t) = e^{t}, 0<t<b
Y(0) = 1, y'(0) = 1, y''(0) = 1
y(t) = e^{t}, b = 1
Problem 3(Homogenous non linear):
Problem 4(Non homogenous linear):
Y'''(t) = 6y^{4}, 0<t<b
Y(0) = 1, y'(0) = 1, y''(0) = 2
PROBLEMS OF DDEs
Problem 5 (Homogenous non linear):
Y(0) = 0, y'(0) = 1, y'(0) = 1
y(t) = e^{t}, b = 1
Problem 6(Homogenous linear):
Y(0) = 1, y'(0) = 1, y'(0) = 1
y(t) = e^{t}, b = 1
Problem 7(Non homogenous non linear)
Problem 8(Non homogenous linear)
Y(0) = 0, y'(0) = 1, y'(0) = 1
y(t) = ln(1+t), b = π
DISCUSSION AND CONCLUSION
In this study, we have derived the RKD method of twostage, third order. This method has been adapted with delay differential equations. The zerostability of the method is proven. The stability polynomial of the method when applied to linear third order DDE is also given. We used the method for solving special third order ODEs as well as DDEs. For DDEs cubic interpolation is used to evaluate the delay terms. Numerical results show that the RKD method is more accurate and requires less computational time compared to the existing RK methods.