**1. Introduction**

The Initial Value Problem (IVP) is

$$y' = f(\mathbf{x}, y), \\ y(\mathbf{x}\_0) = y\_0 \tag{1}$$

with *<sup>x</sup>*<sup>0</sup> <sup>∈</sup> <sup>R</sup> and the vectors *<sup>y</sup>*, *<sup>y</sup>* <sup>∈</sup> <sup>R</sup>*m*. The function *<sup>f</sup>* is defined as *<sup>f</sup>* : <sup>R</sup> <sup>×</sup> <sup>R</sup>*<sup>m</sup>* <sup>→</sup> <sup>R</sup>*m*. Runge–Kutta (RK) pairs are perhaps the most used numerical methods for addressing

(1). They usually presented in a so-called Butcher tableau [1,2] as given below.


In this type of tableau, we have *bT*, ˆ *<sup>b</sup>T*, *<sup>c</sup>* <sup>∈</sup> <sup>R</sup>*<sup>s</sup>* while *<sup>A</sup>* <sup>∈</sup> <sup>R</sup>*s*×*<sup>s</sup>* . Then, the method shares *s* stages and in case that *A* is strictly lower triangular, it is evaluated explicitly. The numerical approximations of the solution step from (*xn*, *yn*) <sup>∈</sup> <sup>R</sup>1+*<sup>m</sup>* to *xn*+<sup>1</sup> <sup>=</sup> *xn* <sup>+</sup> *hn*

by producing two numerical estimations for *<sup>y</sup>*(*xn*+1). Namely, *yn*+<sup>1</sup> <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* and *<sup>y</sup>*ˆ*n*+<sup>1</sup> <sup>∈</sup> <sup>R</sup>*m*, given by

> *s* ∑ *i*=1 *bi fi*,

*i*=1

and

*y*ˆ*n*+<sup>1</sup> = *yn* + *hn s* ∑ ˆ *bi fi*,

*yn*+<sup>1</sup> = *yn* + *hn*

with

$$f\_i = f(\mathfrak{x}\_n + c\_i h\_{n\prime} \mathfrak{y}\_n + h\_n \sum\_{j=1}^{i-1} a\_{ij} f\_j)\_{\prime i}$$

for *<sup>i</sup>* <sup>=</sup> 1, 2, ··· ,*s*. These two approximations *yn*+<sup>1</sup> <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* and *<sup>y</sup>*ˆ*n*+<sup>1</sup> <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* are of algebraic orders *p* and *q* < *p*, respectively. This means that when expanding them in Taylor series, they attain orders *O*(*hp*) and *O*(*hq*), respectively, with *h* being the proper step–length. Then, a local error estimation

$$\varepsilon\_n = h\_n^{p-q-1} \cdot ||y\_{n+1} - \hat{y}\_{n+1}||\_{\infty}.$$

is formed in every step and is combined in an algorithm for changing the steplength.

$$h\_{n+1} = 0.8 \cdot h\_n \cdot \left(\frac{t}{\varepsilon\_n}\right)^{1/p} \cdot t$$

with *t* a small positive number which is set by the user and is named tolerance. The safety factor 0.8 is in common use in such formulas and offers increased reliability to the results. Whenever  *<sup>n</sup>* < *t*, we use the above formula for defining the length of the next step forward. In reverse, when  *<sup>n</sup>* ≥ *t* we again use it without advancing the solution in this case and using instead the value *hn*+<sup>1</sup> as a new trial step *hn*. Information in the issue can be found with details in [3]. As an abbreviation these methods are named RKp(q) airs.

Carl David Tolmé Runge [4] and Martin Wilhelm Kutta [5] introduced the methods bearing their names almost in the turning of the 19th century. For almost 60 years the these methods were implemented with constant step sizes. Richardson extrapolation appeared in the meantime [6] and was used in a kind of step control through doubling and halving [7]. Runge–Kutta pairs appeared 60 years ago. The first series of well-known Runge–Kutta pairs of orders 5(4), 6(5), and 8(7) were presented by Fehlberg [8,9]. In the early 1980s, Dormand and Prince gave their celebrated pairs [10,11].

Non-stiff problems having the form (1) are well suited for being efficiently solved by Runge–Kutta pairs. There is a number of different pairs sharing various orders. We may explain this by the accuracy on demand. Thus, the lesser the accuracy required, the lesser order pairs are chosen. Otherwise, for stringent accuracies at quadruple precision, a high order pair has to be preferred.

Here, we concentrate on RK5(4) pairs which are the first choice when middle tolerances are used. Our special interest is in problems of the form (1) with periodic/oscillatory solutions. In the following, we focus on producing a RK5(4) pair that best address the latter type of problems.

The paper is organized in sections as follows.

