\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
Ninth MSU-UAB Conference on Differential Equations and Computational
Simulations.
\emph{Electronic Journal of Differential Equations},
Conference 20 (2013),  pp. 119--132.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}\setcounter{page}{119}
\title[\hfilneg EJDE-2013/Conf/20/ \hfil Solving oscillatory problems]
{Solving oscillatory problems using a block hybrid
trigonometrically fitted method with two off-step points}

\author[F. F. Ngwane,  S. N. Jator \hfil EJDE-2013/conf/20 \hfilneg]
{Fidele F. Ngwane,  Samuel N. Jator}  % in alphabetical order

\address{Fidele F. Ngwane \newline
Department of Mathematics,
University of south Caroline, Salkehatchie, SC 29488, USA}
\email{fifonge@yahoo.com}

\address{Samuel N. Jator \newline
Department of Mathematics and Statistics,
Austin Peay State University,
Clarksville, TN 37044, USA}
\email{Jators@apsu.edu}

\thanks{Published October 31, 2013.}
\subjclass[2000]{65L05, 65L06}
\keywords{First order system; hybrid method; 
trigonometrically fitted method; \hfill\break\indent 
off-step point; block method}


\begin{abstract}
 A continuous hybrid method using trigonometric basis (CHMTB) with
 two `off-step' points is developed and used to produce three discrete
 hybrid methods which are simultaneously applied as numerical integrators
 by assembling them into a block hybrid trigonometrically fitted method (BHTFM)
 for solving oscillatory initial value problems (IVPs). The stability
 property of the BHTFM is discussed and the performance of the method is
 demonstrated on some numerical examples to show accuracy and efficiency
 advantages.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

In this article, we consider the  first-order differential
equation
\begin{equation} \label{eq-1} 
 y' = f( x,y), \quad y(a) = y_{0}, \quad x \in [a,b],
\end{equation}
 with oscillating solutions where $f:\mathbb{R} \times
\mathbb{R}^m\to \mathbb{R}^m$, $y, y_{0} \in \mathbb{R}^m$.
Oscillatory IVPs frequently arise in areas such as classical
mechanics, celestial mechanics, quantum mechanics, and biological
sciences. Several numerical methods based on the use of polynomial
basis functions have been developed for solving this class of
important problems (see Lambert \cite{Lam1,Lam2},  Hairer et al in
\cite{HWS}, Hairer \cite{HWS10}, and Sommeijer \cite{SJ}). Other methods based
on exponential fitting techniques which take advantage of the special
properties of the solution that may be known in advance have been
proposed (see Simos \cite{SIMOS}, Vanden et al \cite{VAN},  
Vigo-Aguiar et al \cite{RAMOS}, Franco \cite{FRANCO},
Fang et al \cite{FANG}, Nguyen et al  \cite{CONG}, and Jator et al  \cite{JSF}). 
The motivation governing the exponentially-fitted methods is inherent in
the fact that if the frequency or a reasonable estimate of it is
known in advance, these methods will be more advantageous than the
polynomial based methods.

The goal of this article is to construct a CHMTB  which provides
three discrete methods that are combined and applied as a BHTFM which
takes the frequency of the solution as a priori knowledge. The
coefficients of the BHTFM are functions of the frequency and the
stepsize, hence the solutions provided by the proposed method are
exact if \eqref{eq-1} has periodic solutions with known frequencies. We are
motivated to use the hybrid method in order to increase the order of
the method, while preserving good stability properties. Hybrid
methods were first proposed to overcome the Dahlquist \cite{DQ}
barrier theorem whereby the conventional linear multistep method was
modified by incorporating off-step points in the derivation process
(see Gear \cite{GE1}, Gragg et al in  \cite{GS}, Butcher
\cite{BU}, Gupta \cite{GUP}, Lambert \cite{Lam2}, and Kohfeld et al in
\cite{KT1}). These methods were shown to enjoy both higher
order and good stability properties, but included additional
off-grid functions. Gupta \cite{GUP} noted that the design of
algorithms for hybrid methods is more tedious due to the occurrence
of off-step functions which increase the number of predictors needed
to implement the methods. In order to avoid this deficiency, we
adopt a different approach based on a block-by-block implementation
instead of the traditional step-by-step implementation, generally
performed in a predictor-corrector mode.

Hence, we adopted the approach given in Jator et al in \cite{JSF},
where the CHMTB is used to generate a main and two additional
methods which are combined and used as a BHTFM to simultaneously
produce approximations
\[
\{y_{n+\mu}, y_{n+v}, y_{n+1}\}\text{  at a block of points }
\{x_{n+\mu}, x_{n+v}, x_{n+1}\},
\]
$h=x_{n+1}-x_n$, $n=0, \dots, N-1$,  on a partition
$[a,b]$, where $\mu, v \in (0,1)$, $a, b \in  {\mathbb R}$, $h$
is the constant stepsize, $n$ is a grid index and $N>0$ is the
number of steps. Block methods have also been considered by Shampine
and Watts \cite{SHW}. We emphasize that the BHTFM simultaneously
generates approximations $\{y_{n+\mu},  y_{n+v}, y_{n+1}\}$ to the exact
solutions $\{y(x_{n+\mu}), y(x_{n+v}), y(x_{n+1})\}$.

To apply the block method at the next block to obtain
$y_{n+2}$, the only necessary starting value is $y_{n+1}$, and the
loss of accuracy in $y_{n+1}$ does not affect subsequent points,
thus the order of the algorithm is maintained. It is unnecessary to
make a function evaluation at the initial part of the new block.
Thus, at all blocks except the first, the first function evaluation
is already available from the previous block.

The organization of this article is as follows. 
In Section 2, we obtain a
trigonometric basis  representation $U(x)$ for the exact solution
$y(x)$ which is used to generate two discrete methods which are
combined into a BHTFM for solving  \eqref{eq-1}. The analysis and
implementation of the BHTFM are discussed in Section 3. Numerical
examples are given in Section 4 to show the accuracy and efficiency
of the BHTFM. Finally,  we  give some concluding remarks in Section 5.

\section{Development of method}

In this section, our objective is to construct a CHMTB which
produces three discrete methods as by-products. The main method has
the form
\begin{equation} \label{eq-2}
y_{n+1}=y_n + h (\beta_{0}(u) f_n  +  \beta_{1}(u) f_{n+1}
+ \beta_{v}(u) f_{n+v}  +  \beta_{\mu}(u) f_{n+\mu}),
\end{equation}
and the additional methods are given by
\begin{equation}  \label{eq-3}
\begin{gathered}
y_{n+v}=y_n + h (\hat{\beta}_{0}(u) f_n  +  \hat{\beta}_{1}(u)f_{n+1}
+  \hat{\beta}_{v}(u) f_{n+v} +  \hat{\beta}_{\mu}(u) f_{n+\mu}),
\\
y_{n+\mu}=y_n + h (\check{\beta}_{0}(u) f_n  +  \check{\beta}_{1}(u)f_{n+1}
+  \check{\beta}_{v}(u) f_{n+v} +  \check{\beta}_{\mu}(u) f_{n+\mu}),
\end{gathered}
\end{equation}
where $u = w h$, $\beta_{j}(u)$, $\hat{\beta}_{j}(u)$, $\check{\beta}_{j}(u)$,
$\beta_{v}(u)$, $\hat{\beta}_{v}(u)$,  $\check{\beta}_{v}(u)$, $\beta_{\mu}(u)$,
 $\hat{\beta}_{\mu}(u)$   and   $\check{\beta}_{\mu}(u)$, $j=0,1$, are
coefficients that depend on the stepsize and frequency. We note that
$y_{n+v}$  and  $y_{n+\mu}$  are  the numerical approximation to the analytical
solutions $y(x_{n+v})$, and $y(x_{n+\mu})$ respectively and
$f_{n+v}=f(x_{n+v},y_{n+v})$,  $f_{n+\mu}=f(x_{n+\mu},y_{n+\mu})$,
$f_{n+j}=f(x_{n+j},y_{n+j})$ with $j = 0, 1$.
To obtain equations \eqref{eq-2} and \eqref{eq-3}, we
proceed by seeking to approximate the exact solution $y(x)$ on the
interval $[x_n, x_n+h]$ by the interpolating function $U(x)$ of
the form
\begin{equation}   \label{eq-4}
  U(x)= a_0 +a_1 x + a_2 x^2 + a_3 \sin(w x) + a_4 \cos(w x),
\end{equation}
where $a_0, a_1, a_2,  a_3$  and $a_4$  are coefficients that must be
uniquely determined. We then impose that the interpolating function
equation \eqref{eq-4} coincides with the analytical solution at the
end point $x_n$ to obtain the equation
\begin{equation} \label{eq-5}
  U(x_n)=y_n.
\end{equation}
We also require that the function \eqref{eq-4} satisfy
the differential equation \eqref{eq-1}  at the points $x_{n+\mu},  x_{n+v},
x_{n+j}$,   $j=0,1$ to obtain the following set of three equations:
%
\begin{equation}  \label{eq-6}
 U'(x_{n+\mu})=f_{n+\mu},  \quad U'(x_{n+v})=f_{n+v},  \quad
  U'(x_{n+j})=f_{n+j}, \quad j=0,1.
\end{equation}
Equations \eqref{eq-5} and  \eqref{eq-6} lead to a system of five
equations which is solved by Cramer's rule to obtain
$a_j$, $j = 0, 1, 2, 3, 4$. Our continuous CHMTB is constructed by
substituting the values of $a_j$ into equation \eqref{eq-4}.
After some algebraic
manipulation, the CHMTB is expressed in the form
\begin{equation} \label{eq-7}
U(x)=y_n + h(\beta_{0}(w, x) f_n  +  \beta_{1}(w, x) f_{n+1}
+  \beta_{v}(w, x) f_{n+v} + \beta_{\mu}(w, x) f_{n+\mu}),
\end{equation}
 where $w$ is the frequency, $\beta_0(w,x)$,
$\beta_1(w,x)$, $\beta_v(w,x)$,  and $\beta_\mu(w,x)$, are continuous
coefficients. The
continuous method \eqref{eq-7} is used to generate the main method
of the form \eqref{eq-2} and two additional methods of the form
\eqref{eq-3} by choosing, $v=\frac{1}{2}$ and  $\mu=\frac{1}{4}$.
Thus, evaluating \eqref{eq-7} at
$x=\{x_{n+1}, x_{n+\frac{1}{2}},  x_{n+\frac{1}{4}}\}$  and letting
$u = w h$, we obtain the coefficients of \eqref{eq-2} and \eqref{eq-3} as
follows:
\begin{equation}   \label{eq-8}
 \begin{gathered}
\beta_0 = \frac{ \cos(\frac{u}{8}) (\csc(\frac{u}{4})^3)
 \sin(\frac{u}{8}) (u - 2 \sin(\frac{u}{2})) }{2u} \\
\beta_1 = -( \frac{ \cos(\frac{u}{8}) \csc(\frac{u}{4})^3
 \sin(\frac{u}{8}) (-u + 2 \sin(\frac{u}{2}))}{2u})\\
\beta_v = -( \frac{\cos(\frac{u}{8}) \csc(\frac{u}{4})^3
 \sin(\frac{u}{8}) (u \cos(\frac{u}{2}) - 2 \sin(\frac{u}{2}))}{u}),
\end{gathered}
\end{equation}
and 
\begin{equation}  \label{eq-9}
\begin{gathered}
\hat{\beta}_0 = \frac{ (\csc(\frac{u}{8})^2) (u - 4 \sin(\frac{u}{4})) }{8 u}
\\
\hat{\beta}_v  = \frac{(\csc(\frac{u}{8})^2) (u - 4 \sin(\frac{u}{4}))} {8 u}
\\
\hat{\beta}_\mu  = -(\frac{ (\csc(\frac{u}{8})^2) (u \cos(\frac{u}{4})
- 4 \sin(\frac{u}{4}))}{4 u})
\\
\check{\beta}_0 = \frac{(\csc(\frac{u}{4})^3)
\sin(\frac{u}{8}) (8 u \cos(\frac{u}{8}) + 3 u \cos(\frac{(3 u)}{8})
 -  8 (2 \sin(\frac{(3 u)}{8}) + \sin(\frac{(5 u)}{8})))} {16 u}
\\
\check{\beta}_1  =  -(\frac{(\csc(\frac{u}{4})^3) (u \cos(\frac{u}{8})
- 8 \sin(\frac{u}{8})) \sin(\frac{u}{8})} {16 u})
\\
\check{\beta}_v  = \frac{(3 + 3 \cos(\frac{u}{4}) + \cos(\frac{u}{2}))
 (  \csc(\frac{u}{4})^3) (u \cos(\frac{u}{8}) - 8 \sin(\frac{u}{8}))
 \sin(\frac{u}{8})} {8 u}
\\
\check{\beta}_\mu  = -( \frac{(\cos(\frac{u}{8})^2)
(\csc(\frac{u}{4})^3)  \sin(\frac{u}{8}) (3 u \cos(\frac{u}{8})
+ 3 u \cos(\frac{(3 u)}{8}) - 16 \sin(\frac{(3 u)}{8}))} { 4 u}),
\end{gathered}
\end{equation}
with $\hat{\beta}_1  = 0 $  and $ \beta_\mu = 0$.

\section{Error analysis and stability}

\subsection{Local truncation error} 
The Taylor  series is used for  small values of $u$ (see
Simos \cite{SIMOS}). Thus the coefficients  in equation \eqref{eq-8}
can be expressed as
\begin{gather*} %10
\beta_0 = \frac{1}{6} + \frac{u^2}{720} + \frac{u^4}{80640} 
+ \frac{u^6}{9676800} + \frac{u^8}{1226244096}
 + \frac{691 u^{10}}{111588212736000}+ \dots
\\
\beta_1 = \frac{1}{6} + \frac{u^2}{720} + \frac{u^4}{80640} 
+ \frac{u^6}{9676800} + \frac{u^8}{1226244096} 
+ \frac{691u^{10}}{111588212736000} + \dots
\\
\beta_v = \frac{2}{3} - \frac{u^2}{360} - \frac{u^4}{40320} 
- \frac{u^6}{4838400} - \frac{u^8}{613122048} - \frac{691u^{10}}{55794106368000}
  +  \dots
\end{gather*} 
and  the coefficients in equation \eqref{eq-9} can be expressed as
\begin{gather*}   %11
\begin{aligned}
\hat{\beta}_0 &= \frac{1}{12} + \frac{u^2}{5760} + \frac{u^4}{2580480} 
+ \frac{u^6}{1238630400} + \frac{u^8}{627836977152} \\
&\quad + \frac{691 u^{10}}{228532659683328000}   +  \dots
\end{aligned}\\
\begin{aligned}
\hat{\beta}_v &= \frac{1}{12} + \frac{u^2}{5760} + \frac{u^4}{2580480} 
+ \frac{u^6}{1238630400} + \frac{u^8}{627836977152} \\
&\quad + \frac{691 u^{10}}{228532659683328000}   +  \dots
\end{aligned}\\
\begin{aligned}
\hat{\beta}_\mu &= \frac{1}{3} - \frac{u^2}{2880} - \frac{u^4}{1290240} 
- \frac{u^6}{619315200} - \frac{u^8}{313918488576} \\
&\quad - \frac{691 u^{10}}{114266329841664000}   +  \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_0  &=  \frac{37}{384} + \frac{67u^2}{184320} 
+ \frac{401u^4}{165150720} + \frac{1649u^6}{79272345600} 
+ \frac{1711u^8}{9132174213120} \\
&\quad + \frac{12094087 u^{10}}{7313045109866496000}   
+  \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_1  &=  \frac{1}{384} + \frac{13 u^2}{184320} 
+ \frac{37 u^4}{33030144} + \frac{1091 u^6}{79272345600} 
+ \frac{2087 u^8}{14350559477760} \\
&\quad + \frac{10191073 u^{10}}{7313045109866496000} 
 +  \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_v & = -\frac{7}{192} + \frac{7 u^2}{46080} 
- \frac{11 u^4}{11796480} - \frac{29 u^6}{1415577600} 
- \frac{12503 u^8}{50226958172160}\\
&\quad  - \frac{659969 u^{10}}{261180182495232000}  
  +  \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_\mu & =  \frac{3}{16} - \frac{3 u^2}{5120} - \frac{3 u^4}{1146880} 
- \frac{31 u^6}{2202009600} - \frac{13 u^8}{155021475840} \\
&\quad - \frac{11747 u^{10}}{22571126882304000}  +  \dots.
\end{aligned}
\end{gather*} 

Thus, the Local Truncation Errors (LTEs) for methods
\eqref{eq-2} and \eqref{eq-3} are given by
\begin{equation} %12
\begin{gathered}
\text{LTE}\eqref{eq-2} = -\frac{h^5}{2880}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)),
\\
\text{LTE}\eqref{eq-3}  = -\frac{h^5}{92160}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)),
\\
\text{LTE}\eqref{eq-3}  = -\frac{53h^5}{1474560}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)).
\end{gathered}
\end{equation}


\begin{remark} \label{rmk3.1} \rm 
We noticed that as $u \to 0$, the method (2)
reduces to the fourth-order method of  Gragg and Stetter \cite{GS},
which is based on polynomial basis functions.
\end{remark}

\subsection{Stability}

The BHTFM is constructed by combining equations \eqref{eq-2} and
\eqref{eq-3}, where the coefficients of the method are explicitly
given by equations \eqref{eq-8} and \eqref{eq-9}. We then define the
block-by-block method for computing vectors $Y_0, Y_1, \cdots $ in
sequence (see \cite{FT}). Let the $\eta-$vector ($\eta = 3$ is the
number of points within the block) $Y_\gamma, Y_{\gamma -1},  F_\gamma,
  F_{\gamma - 1}$ be given as 
\begin{gather*}
Y_\gamma = (y_{n+ \frac{1}{4}}, y_{n+ \frac{1}{2}}, y_{n+1})^T, \quad
Y_{\gamma-1} = (y_{n - \frac{1}{4}}, y_{n - \frac{1}{2}}, y_n)^T, \\
F_\gamma = ( f_{n+ \frac{1}{4}},  f_{n+ \frac{1}{2}}, f_{n+1})^T, \quad
F_{\gamma-1} = (f_{n - \frac{1}{4}},  f_{n - \frac{1}{2}},  f_n)^T,
\end{gather*}
then the $1$-block $3$-point method for equation \eqref{eq-1} is
given by
\begin{equation}\label{eq-13}
Y_\gamma = \sum_{i=1}^1 A^{(i)}Y_{\gamma-i} +  \sum_{i=0}^1
B^{(i)}F_{\gamma-i},
\end{equation}
where $A^{(i)}, B^{(i)}$, $i = 0, 1$ are  $3 \times 3 $ matrices
whose entries are given by the coefficients of \eqref{eq-2} and \eqref{eq-3}.

\subsection*{Zero-stability} 
\begin{definition}\label{df-1} \rm
The block method \eqref{eq-13} is zero stable provided the roots $R_j$, 
$j=1, 2, 3$ of the  first characteristic polynomial $\rho(R)$ specified
by
\begin{equation}\label{eq-14}
\rho(R) =\det  \Big[  \sum_{i=0}^1 A^{(i)}R^{1-i} \Big]  =
0, \quad A^{(0)} = -I,
\end{equation}
satisfies $|R_j| \leq 1$, $j= 1, 2, 3$   and  for those roots with $|R_j| =1$, 
the multiplicity does not exceed $1$ (see \cite{FT}).
\end{definition}

\subsection*{Consistency of BHTFM} 
We note that the block method \eqref{eq-13} is consistent as  it has 
order $p>1$. We see from equation \eqref{eq-14} and  definition \eqref{df-1} 
that the block method \eqref{eq-13} is zero-stable since $\rho(R)=R^2(R-1)=0$
satisfies $|R_j| \leq 1$, $j=1, 2, 3$, and for those roots with 
$|R_j| = 1$, the multiplicity does not exceed $1$. The block method
\eqref{eq-13}  is thus convergent, as zero-stability plus consistency
equals  convergence.

\subsection*{Linear stability of the BHTFM}   
To analyze the linear stability of the BHTFM, we  apply the method to the test
equation $y' =\lambda y$, where $\lambda$ is expected to run
through the (negative) eigenvalues of the Jacobian matrix. Then  an
application of \eqref{eq-13} to the test equation yields
\begin{equation}\label{eq-15}
Y_\gamma = M(q;u)Y_{\gamma-i},  \quad q = h \lambda,  \quad u = w h,
\end{equation}
where
\[
M(q;u) = (A^{(0)}- qB^{(0)})^{-1}   (A^{(1)}+ qB^{(1)}),
\]
where the matrix $M(q;u)$ is the  amplification matrix which determines
 the stability of the method.

\begin{definition}\label{df-2} \rm
A region of stability is a region  in the $q-u$ plane, throughout  
which $|\rho(q;u)| \leq 1$, where $\rho(q;u)$ is the spectral radius of
 $M(q;u)$.
\end{definition}

It is easily seen that the  eigenvalues of  $M(q;u)$ are
$\lambda_1 =0$,  $\lambda_2 =0$, and 
\begin{align*}
\lambda_3 &=
 \Big(-16 (2 + q) (q^2 + u^2) \cos(\frac{u}{8}) + (6q^3 + 16u^2 + 6qu^2 
+ q^2(16 + u^2)) \cos(\frac{3u}{8}) \\
&\quad +   16 q^2 \cos(\frac{5u}{8}) 
 + 10q^3 \cos(\frac{5u}{8}) + 16u^2 \cos(\frac{5u}{8}) 
 +  10qu^2 \cos(\frac{5u}{8}) \\
&\quad + 3q^2u^2\cos(\frac{5u}{8})
  +   q^3u \sin(\frac{3u}{8}) +   3q^3u \sin(\frac{5 u}{8})\Big)\\
&\quad\div   
   \Big(16 (-2 + q) (q^2 + u^2) \cos(\frac{u}{8}) + (-10 q^3 
 + 16 u^2 - 10 q u^2 \\
&\quad + q^2 (16 + 3 u^2)) \cos(\frac{3 u}{8}) + 16 q^2 \cos(\frac{5 u}{8})  
 - 6 q^3 \cos(\frac{5 u}{8}) +   16 u^2 \cos(\frac{5 u}{8}) \\
&\quad - 6 q u^2 \cos(\frac{5 u}{8}) 
 + q^2 u^2 \cos(\frac{5 u}{8}) 
 - 3 q^3 u \sin(\frac{3 u}{8}) - q^3 u \sin(\frac{5 u}{8})\Big).
 \end{align*}
 
\begin{figure}[ht]
  \begin{center}
  \includegraphics[width=0.4\textwidth]{fig1} % StaRegP3
\end{center}
  \caption{The shaded region represents the truncated region of absolute 
stability}
  \label{StaReg}
\end{figure}
%
\begin{figure}[ht]
  \begin{center}
     \includegraphics[width=0.4\textwidth]{fig2} % ZeroPolesP3 
\end{center}
  \caption{$\lambda_{3}$ has three zeros$(\square)$ and  no poles$(+)$ in $\mathbb{C}^-$, hence the BHTFM is $A_0$-stable.}
  \label{ZerosPoles}
\end{figure}


We observed that in the $q-u$ plane the BHTFM  is stable for 
$ q \in [-100, 100]$, and $u \in [-\pi, \pi ]$.
Figure \ref{StaReg} is a plot of the stability region. Figure \ref{ZerosPoles} 
shows the zeros and poles of $\lambda_3$.


\subsection{Implementation} 
We emphasize that methods \eqref{eq-2} and \eqref{eq-3} are combined 
to form the block method \eqref{eq-13}, which is implemented to solve 
\eqref{eq-1} without requiring
starting values and predictors. For instance, if we let $n=0$ and
$\gamma = 1$ in \eqref{eq-13}, then 
$(y_{\frac{1}{4}}, y_{\frac{1}{2}}, y_{1})^{T}$ are
simultaneously obtained over the sub-interval $[x_{0},x_{1}]$, as
$y_{0}$ is known from the IVP. Similarly, if $n = 1$ and $\gamma = 2$, then
$(y_{\frac{5}{4}}, y_{\frac{3}{2}}, y_{2} )^{T}$ are simultaneously obtained 
over the sub-interval $[x_{1},x_{2}]$, as $y_{1}$ is known from the previous
block, and so on; until we reach the final sub-interval
$[x_{N-1},x_{N}]$. We note that for linear problems, we solve
\eqref{eq-1} directly using the feature solve[~~] in Matlab, while
nonlinear problems were solved by implementing the Newton's method
in Matlab enhanced by the feature fsolve[~~].


\section{Numerical examples}

In this section, we give numerical examples to illustrate the
accuracy (small errors) and efficiency (fewer number of function
evaluations (FNCs)) of the BHTFM. We find the approximate solution
on the partition $\pi_{N}$, where 
\begin{equation}
\label{pi} 
\pi_{N}:a=x_{0} <x_{1}
< x_{2} < \dots < x_n < x_{n+1} < \dots < x_{N}=b, 
\end{equation}
and we give the
errors at the endpoints calculated as Error=$y_{N}-y(x_{N})$. We
note that the method requires only three function evaluations per step
and in general requires $(3N + 1)$ FNCs on the entire interval.  All
computations were carried out using a written code in Matlab.


\begin{example}\label{Example-1} \rm
We consider the  following inhomogeneous IVP by Simos  \cite{SIMOS}.
\[
y'' = - 100y + 99 \sin(x),  \quad  y(0) = 1, \quad y'(0) =11, \quad
 x \in [0, 1000] 
\]
where the analytical solution is given by
 \[
\text{Exact: } y(x) = \cos(10x) + \sin(10x) + \sin(x). 
\]
\end{example}

\begin{table}[ht]
\begin{center} 
\begin{tabular}{ccccccccccccccccc}
\hline
\multicolumn{3}{c}{BHTFM }  & & &  \multicolumn{2}{c}{Simos \cite{SIMOS}}  \\
 N & |Error| &FNCs         & & &  |Error|  &FNCs \\
\hline
1000 & $1.2 \times 10^{-3}$&6002      &  & &   $1.4 \times 10^{-1}$& 8000\\
2000 & $1.2 \times 10^{-3}$&12002      &  & &   $3.5 \times 10^{-2}$& 16000\\ 
4000 & $1.4 \times 10^{-5}$&24002      &  & &   $1.1 \times 10^{-3}$& 32000\\
8000 & $1.5 \times 10^{-7}$&48002      &  & &   $8.4 \times 10^{-5}$& 64000\\
16000 & $8.7 \times 10^{-9}$&96002      &  & &   $5.5 \times 10^{-6}$& 128000\\
32000 & $1.1 \times 10^{-9}$&192002      &  & &   $3.5 \times 10^{-7}$& 256000\\
\hline
\end{tabular} 
\end{center}
\caption{Results, with $\omega = 10$, for example
\eqref{Example-1}} \label{table-1}
\end{table}

The exponentially-fitted method in Simos \cite{SIMOS} is fourth order 
and hence comparable to our method, BHTFM. We see from Table \eqref{table-1} 
that BHTFM is more efficient than the  method in Simos   \cite{SIMOS}.
 We also compare the computational efficiency of the two methods by 
considering the FNCs over $N$ integration steps for each method. 
Our method, BHTFM,   requires only $3N + 1$ function evaluations in $N$ steps 
 compared to $4N$  function evaluations in $N$  steps  for the method 
in Simos  \cite{SIMOS}. Hence for this example, BHTFM  performs better.

\begin{figure}[ht]
  \begin{center}
   \includegraphics[width=0.7\textwidth]{fig3} % EffCurExample1
\end{center}
  \caption{Efficiency curves for  example \eqref{Example-1}}
  \label{EC-1}
\end{figure}



\begin{example}\label{Example-2} \rm
 We consider the IVP (see Vigo-Aguiar et al \cite{RAMOS})
\[
 y'' + K^2y = K^2x,\quad y(0) = 10^{-5},\quad y'(0) = 1 - K 10^{-5} \cot(K),\quad
 x \in [0,100]
\]
where  $ K = 314.16$, and we choose $\omega = 314.16$. The analytical 
solution is given by
\[
\text{Exact: }  y(x) = x + 10^{-5} ( \cos(Kx) - \cot(K) \sin(Kx)).
\]
\end{example}


\begin{table}[ht]
\begin{center} 
\begin{tabular}{cccccccccccccccc}
\hline
\multicolumn{3}{c}{BHTFM }  &  & \multicolumn{3}{c}{CHEBY24}  \\
 N &|Error| &FNCs         & & N&  |Error|  &FNCs \\
\hline
 9 & $5.07 \times 10^{-11}$&48        & &  9 & $1.84 \times 10^{-11}$& 450\\
20 & $9.17 \times 10^{-12}$&122       & &   &   & \\
 40 & $4 \times 10^{-15}$&242        & &   &  & \\
\hline
\end{tabular} 
\end{center}
\caption{Results, with $\omega = 314.16$, for example
\eqref{Example-2}  on $[0,100].$} \label{table-2a}
\end{table}


\begin{table}[ht]
\begin{center} 
\begin{tabular}{cccccccccccccccc}
\hline
\multicolumn{3}{c}{BHTFM }  &  & \multicolumn{3}{c}{CHEBY1}  \\
 N &|Error| &FNCs         & & N&  |Error|  &FNCs \\
\hline
 2 & $4.13 \times 10^{-17}$&14        & &  1 & $1 \times 10^{-16}$& 8\\
\hline
\end{tabular}
\end{center}
\caption{Results, with $\omega = 314.16$, for example
\eqref{Example-2}  on $[0,1]$} \label{table-2b}
\end{table}

This problem demonstrates the performance of BHTFM on  a well-known 
oscillatory problem. We compare the results from BHTFM  with  the
 Dissipative Chebyshev exponential-fitted methods, CHEBY24   and  
CHEBY1  used in Vigo-Aguiar et al \cite{RAMOS}.  We see that BHTFM   
uses fewer number of function evaluations with better accuracy than CHEBY24 
that is designed to use fewer number of steps.
Integrating in the interval $[0,1]$ with a stepsize equal to half the total 
length of the interval, we obtain an error of order $10^{-17}$. 
Hence BHTFM  is a  more efficient integrator.
We note that compared with the methods CHEBY24   and CHEBY1 which use stepsizes 
considerably larger than those used in multistep methods,  BHTFM is very 
competitive to CHEBY1  and superior to CHEBY24.


\begin{example}\label{Example-3} \rm
We consider the  nonlinear perturbed system on the range $[0,10]$, 
 with  $\epsilon = 10^{-3}$  ( see Fang et al \cite{FANG}).
\begin{gather*}
y_1'' + 25y_1 + \epsilon(y_1^2 + y_2^2) 
 = \epsilon \phi_1(x), \quad y_1(0) = 1, \quad y_1'(0) = 0,\\
y_2'' + 25y_2 + \epsilon(y_1^2 + y_2^2) = \epsilon \phi_2(x), \quad
 y_2(0) = \epsilon, \quad y_2'(0) = 5,
\end{gather*}
where
\begin{gather*}
\phi_1(x) = 1 + \epsilon^2 + 2\epsilon \sin(5x+x^2) + 2 \cos(x^2) 
+ (25 - 4x^2) \sin(x^2) ,\\  
\phi_2(x) = 1 + \epsilon^2 + 2\epsilon \sin(5x+x^2) - 2 \sin(x^2) 
+ (25 - 4x^2) \cos(x^2),
\end{gather*}
 and the  exact solution is given by
\[\text{Exact: } y_1(x)= \cos(5x) + \epsilon \sin(x^2),\quad
y_2(x)= \sin(5x) + \epsilon \cos(x^2),
\]
represents a periodic motion of constant frequency with small perturbation 
of variable frequency.
\end{example}

\begin{table}[ht]
\begin{center} 
\begin{tabular}{cccccccccccccccc}
\hline
 \multicolumn{2}{c}{BHTFM} &    \multicolumn{2}{c}{ARKN(5)}&     
\multicolumn{2}{c}{TFARKN(5)}   \\
 $N$ &$-\log_{10}(\rm Err)$ &    $N$(rejected) & $-\log_{10}(\rm Err)$&     
 $N$(rejected)& $-\log_{10}(\rm Err)$  \\
\hline
50 &  $4.04$ &  $42$(15) & $2.82$&  $29$(6)& $2.78$  \\
90 &  $5.04$ &  $86$(7)  & $4.96$&  $88$(9)& $5.33$  \\
170 &  $6.07$ & $260$(5) & $7.16$&  $262$(8)& $7.85$ \\
\hline
\end{tabular} 
\end{center}
\caption{Results, with $\omega = 5$, for example
\eqref{Example-3}} \label{table-3}
\end{table}


We use this problem  to demonstrate the performance of the BHTDA on a nonlinear
perturbed system. This problem was also solved by Fang et al \cite{FANG} 
using a variable stepsize fifth-order trigonometrically fitted
 Runge-Kutta-Nystr\"om method TFARKN5(3)
and a fifth-order Runge-Kutta-Nystr\"om method (ARKN5(3)) which was constructed
by Franco \cite{FRANCO}.  We compare the maximum global error 
$({\rm Err} = \max|y(x)-y|)$
for the three methods in Table \ref{table-3}.  We remark that  TFARKN5(3) 
and ARKN5(3) are expected to perform better because they are exact 
when the solution involves a linear combination of trigonometric 
functions as well as implemented as a variable-step
method. However, BHTFM which is implemented using a fixed step-size
is highly competitive to them.

\begin{figure}[ht]
  \begin{center}
   \includegraphics[width=0.7\textwidth]{fig4} % EffCurExample3
\end{center}
  \caption{Efficiency curves for  example \eqref{Example-3}}
  \label{EC-3}
\end{figure}



\begin{example}\label{Example-4} \rm
We consider the  nonlinear Duffing equation which was  also solved by 
Simos \cite{SIMOS} and Ixaru et al \cite{IXARU}
\[
y'' + y + y^3 = B \cos(\Omega x), \quad y(0) = C_0, \quad y'(0) = 0. 
\]
The analytical solution is given by
\[
\text{Exact:}  y(x) = C_1 \cos(\Omega x)  + C_2 \cos(3\Omega x)  
+  C_3 \cos(5\Omega x)  + C_4 \cos(7\Omega x),
\]
where $\Omega = 1.01$, $B = 0.002$, $C_0 = 0.200426728069$,
$C_1 = 0.200179477536$,  $C_2 = 0.246946143 \times 10^{-3}$,
$C_3 = 0.304016 \times 10^{-6}$, $C_4 = 0.374 \times 10^{-9}$.
  We choose $\omega = 1.01$
\end{example}

\begin{table}[ht]
\begin{center} 
\begin{tabular}{cccccccccccccccc}
\hline
 \multicolumn{2}{c}{BHTFM}   & &    \multicolumn{2}{c}{Simos}&      
 & \multicolumn{2}{c}{Ixaru et al}   \\
 $N$ &|Error| &      & $N$ & |Error|&     & $N$& |Error|   \\
\hline
 $150$ &$ 1.3 \times 10^{-3}$ &      & $300$ & $1.7 \times 10^{-3}$&     
 & $300$& $1.1 \times 10^{-3}$   \\
 $300$ &$ 5.6 \times 10^{-5}$ &      & $600$ & $1.9 \times 10^{-4}$&     
 & $600$& $5.4 \times 10^{-5}$   \\
 $600$ &$ 3.2 \times 10^{-6}$ &      & $1200$ & $1.4 \times 10^{-5}$&    
 & $1200$& $1.9 \times 10^{-6}$   \\
 $1200$ &$ 1.7 \times 10^{-7}$ &      & $2400$ & $8.7 \times 10^{-7}$&  
  & $2400$& $6.2 \times 10^{-8}$  \\
\hline
\end{tabular}
\end{center}
\caption{Results, with $\omega = 1.01$, for example
\eqref{Example-4}} \label{table-4}
\end{table}

We compare  the end-point global errors for BHTFM    with the fourth 
order methods in Simos \cite{SIMOS} and Ixaru et al \cite{IXARU}.
We see from Table \ref{table-4} that the results produced by BHTFM 
are better than those given Simos \cite{SIMOS}, as  BHTFM   
produces  better error magnitude while using only half the number 
of steps needed by Simos \cite{SIMOS}. BHTFM  is very competitive  
to the method used by Ixaru et al \cite{IXARU}.


\begin{figure}[ht]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig5} % EffCurExample4
\end{center}
  \caption{Efficiency curves for  example \eqref{Example-4}}
  \label{EC-4}
\end{figure}


\begin{example}\label{Example-5} \rm
 A nearly sinusoidal problem.
We consider the following IVP on the range $0\leq t \leq 10$, 
(see  Nguyen et al  \cite[p. 205]{CONG})
\begin{gather*}
y_1'= -2y_1 + y_2 + \sin(t),\quad  y_1(0) = 2, \\\
y_2'= -(\beta +2)y_1 + (\beta+1)y_2 + \sin(t) - \cos(t),\quad y_2(0) = 3.
\end{gather*}
We choose $\beta = -3$ and $\beta =-1000$  to illustrate the phenomenon 
of stiffness. Given the initial conditions $y_1(0) = 2$ and $y_2(0) = 3$, 
the exact solution is $\beta$-independent and is given by
\[
\text{Exact: } y_1(t)=  2 \exp(-t) + \sin(t),\quad
 y_2(t)=  2 \exp(-t) + \cos(t).
\]
\end{example}

\begin{table}[ht]
\begin{center} 
\begin{tabular}{cccccccccccccccc}
\hline
\multicolumn{3}{c}{BHTFM with $(\beta = -3)$} &  
& \multicolumn{3}{c}{ Nguyen et al \cite{CONG} with $(\beta = -3)$}  \\
 N &|Error| &FNCs         & &  N & |Error|  &FNCs \\
\hline
 6 & $8.9 \times 10^{-6}$&38      &  &   10 & $5.4 \times 10^{-6}$& 47\\
 10 & $9.0 \times 10^{-7}$&62      &  &   19 & $8.3 \times 10^{-8}$& 88\\
 19 & $5.8 \times 10^{-8}$&116      &  &   23 & $4.5 \times 10^{-4}$& 113\\
\hline
\end{tabular} 
\end{center}
\caption{Results, with $\omega = 1$, for example
\eqref{Example-5} with $\beta =-3.$} \label{table-5a}
\end{table}


\begin{table}[ht]
\begin{center}
\begin{tabular}{cccccccccccccccc}
\hline
\multicolumn{3}{c}{BHTFM with $(\beta = -1000)$}  &  & 
\multicolumn{3}{c}{Nguyen et al \cite{CONG} with $(\beta = -1000)$}  \\
 N &|Error| &FNCs    & & N& |Error|  &FNCs \\
\hline
 6 & $8.9 \times 10^{-6}$&38      &  &  13 & $1.0 \times 10^{-6}$& 61\\
 10 & $9 \times 10^{-7}$&62      &  &  16 & $1.6 \times 10^{-7}$& 76\\
 13 & $2.9 \times 10^{-7}$&80      &  &  21 & $7.0 \times 10^{-8}$& 98\\
 16 & $1.1 \times 10^{-7}$&99      &  &   &    \\
 21 & $3.8 \times 10^{-8}$&128      &  & &   &    \\
\hline
\end{tabular} 
\end{center}
\caption{Results, with $\omega = 1$, for example
\eqref{Example-5} with $\beta = -1000$.} \label{table-5b}
\end{table}

This  example is chosen to demonstrate the performance of BHTFM  on 
 stiff problems. We compute the  solutions to Example \eqref{Example-5}   
with $\beta =-3,  -1000$.  We obtain comparable or better absolute errors 
than Nguyen et al (\cite{CONG}). This efficiency is achieved  using  
fewer number of blocks and less function evaluations than  Nguyen et al  
(\cite{CONG}). For example when $\beta =-3$, our method generates a solution 
with error magnitude $10^{-6}$ involving just $6$ blocks  and $38$ function 
evaluations,  whereas \cite{CONG}  attains the same error  magnitude  
using $10$ blocks  and $47$ function evaluations.
When $\beta =-1000$ and using $10$ blocks, BHTFM generates a solution 
with error magnitude $10^{-7}$ involving $62$ function evaluations, 
 whereas \cite{CONG}  attains an error  magnitude of $10^{-6}$ while  
using $13$  blocks.  We see that BHTFM performs better using  fewer blocks  
and is competitive to Nguyen et al (\cite{CONG}) which is of order six and 
is thus expected to do better.

\begin{figure}[ht]
  \begin{center}
\includegraphics[width=0.7\textwidth]{fig6} % EffCurExample5a
\end{center}
  \caption{Efficiency curves for  example \eqref{Example-5} with $\beta =-3$}
  \label{EC-5a}
\end{figure}

\begin{figure}[ht]
  \begin{center}
    \includegraphics[width=0.7\textwidth]{fig7} % EffCurExample5b
\end{center}
  \caption{Efficiency curves for  example \eqref{Example-5} 
with $\beta = -1000$}
  \label{EC-5b}
\end{figure}

\begin{example}\label{Example-6} \rm  Linear Kramarz problem.
We consider the following second-order IVP, (see  Nguyen et al 
\cite[p. 204]{CONG})
\begin{gather*}
y''(t)= \begin{pmatrix} 2498 & 4998 \\
-2499  & -4999 \end{pmatrix} y(t),  \quad
y(0)= \begin{pmatrix}   2\\
-1 \end{pmatrix},  \quad
y'(0) = \begin{pmatrix}   0 \\
 0 \end{pmatrix}, \quad 0 \leq t \leq 100.
\end{gather*} 
Exact: $y(t)= ( 2 \cos(t), -\cos(t))^T$.
\end{example}

We use this example to show the efficiency of BHTFM on linear systems. 
Nguyen et al \cite{CONG} used the "trigonometric implicit
Runke-Kutta", TIRK3,  method  to solve the above linear Kramarz problem.  
Clearly, BHTFM performs better  as seen in Table \ref{table-6}.

\begin{table}[ht]
\begin{center}
\begin{tabular}{cccccccccccccccc}
\hline
\multicolumn{3}{c}{BHTFM} & & \multicolumn{3}{c}{Nguyen et al \cite{CONG}}  \\
 N &|Error| &FNCs         &  &N& |Error|  &FNCs \\
\hline
10&  $8.3\times 10^{-15}$&144   &  &     73&$3.3\times 10^{-12}$ & 327  \\
30&  $5 \times 10^{-14}$&364    &  &      142 &$0.9\times 10^{-11}$ & 707  \\
40&  $7.2 \times 10^{-14}$&484   &  &        170 &$3.7\times 10^{-12}$ & 811\\
43&  $9.5\times 10^{-14}$&520    &  &       -& -&  -& \\
\hline
\end{tabular} 
\end{center}
\caption{Results, with $\omega = 1$, for example
\eqref{Example-6}} \label{table-6}
\end{table}

\begin{figure}[ht]
  \begin{center}
 \includegraphics[width=0.7\textwidth]{fig8} % EffCurExample6
\end{center}
  \caption{Efficiency curves for  example \eqref{Example-6}}
  \label{EC-6}
\end{figure}


\subsection{Estimating the frequency}\label{FreEst}
A preliminary testing indicates that a good estimate of the
frequency can be obtained by demanding that $LTE(2)-LTE(3)=0$, and
solving for the frequency quadratically. That is, solve for $\omega$ given
that
\[
\Big(-\frac{h^5}{2880}(w^2y^{(3)}(x_n) + y^{(5)}(x_n))\Big) - 
 \Big(-\frac{h^5}{92160}(w^2y^{(3)}(x_n) + y^{(5)}(x_n))\Big) = 0,
\]
 where $y^{(j)}$, $j = 2, \dots, 5$ are derivatives that can
be obtained analytically from the given differential equation or
could be calculated via the Taylor series expansion. We used this
procedure to calculate $\omega$ for the problem given in example 
\eqref{Example-1} and obtained
\[
\omega =\pm   \sqrt(9091/91)\approx \pm 9.99505,
\]
which  approximately gives  the known frequency $ \omega = 10.0$. Hence,
this procedure is interesting and will be seriously considered in
our future research.

\subsection*{Conclusion}
We have proposed a BHTFM  for solving periodic IVPs. The BHTFM is
$A_0$-stable and hence, an excellent candidate for solving stiff IVPs.
This method  has only two 'off-step' points and has the advantages of
being self-starting, having good accuracy with order $4$, and
requiring only three functions evaluation  at each integration step.
We have presented representative numerical examples that are linear,
non-linear, stiff and highly oscillatory. These examples show that the
BHTFM is more accurate and  efficient than those in Nguyen et al
\cite{CONG}, Simos \cite{SIMOS}, Vigo-Aguiar et al \cite{RAMOS},   
and Fang \cite{FANG}. Details of
the numerical results are displayed in Tables \ref{table-1}, 
\ref{table-2a}, \ref{table-2b}, \ref{table-3}, \ref{table-4},
\ref{table-5a}, \ref{table-5b}, \ref{table-6}  and the
efficiency curves are presented in Figures \ref{EC-1}, \ref{EC-3},  
\ref{EC-4}, \ref{EC-5a}, \ref{EC-5b}, \ref{EC-6} . Our
future research will incorporate a technique for accurately
estimating the frequency as suggested in subsection \ref{FreEst} as well as
implementing the method in a variable step mode.

\begin{thebibliography}{00}

\bibitem{BU} J. C. Butcher; 
\emph{A modified multistep method for the numerical
integration of ordinary differential equations}, J. Assoc. Comput. Mach. 12 (1965) 124-135.

\bibitem{DQ} G. G. Dahlquist; 
\emph{Numerical integration of ordinary differential equations}, 
Math. Scand. 4 (1956) 69-86.

\bibitem{FANG} Y. Fang, Y. Song and X. Wu; 
\emph{A robust trigonometrically fitted embedded pair for perturbed oscillators},
 J. Comput. Appl. Math., 225 (2009)  347-355.

\bibitem{FT} S. O. Fatunla; 
\emph{Block methods for second order IVPs}, Intern. J. Comput. Math. 41 (1991). 
55 - 63.

\bibitem{FRANCO} J. M.  Franco; 
\emph{Runge-Kutta-Nystr\"{o}m methods adapted to the numerical intergration 
of perturbed oscillators}, Comput. Phys. Comm.  147 (2002)  770-787.

\bibitem{GE1} C. W. Gear; 
\emph{Hybrid methods for initial value problems in ordinary differential 
equations}, SIAM J. Numer. Anal. 2 (1965) 69-86.

\bibitem{GS}W. Gragg and H. J. Stetter; 
\emph{Generalized multistep predictor-corrector methods}, 
J. Assoc. Comput. Mach. 11 (1964) 188-209.

\bibitem{GUP}G. K. Gupta; 
\emph{Implementing second-derivative multistep methods using Nordsieck 
polynomial representation}, Math. Comp. 32 (1978) 13-18.

\bibitem{HWS} E. Hairer and  G. Wanner; 
\emph{Solving Ordinary Differential Equations II}, Springer, New York, 1996.

\bibitem{HWS10} E. Hairer; 
\emph{A One-step Method of Order 10 for $y''=f(x,y)$}, IMA J. Numer. 
Anal. 2, (1982)  83-94.

\bibitem{IXARU} L. Ixaru, and  G. V. Berghe; 
\emph{Exponential fitting},  Kluwer, Dordrecht,  Netherlands (2004).

\bibitem{JSF} S. N. Jator, S. Swindle, and R. French;
 \emph{Trigonometrically fitted block Numerov type method for 
$y''=f(x, y, y')$}, Numerical Algorithms (2012), DOI 10.1007/s11075-012-9562-1.

\bibitem{KT1} J. J. Kohfeld and G. T. Thompson; 
\emph{Multistep methods with modified predictors and correctors}, 
J. Assoc. Comput. Mach. 14 (1967) 155-166.

\bibitem{Lam1} J. D. Lambert; 
\emph{ Numerical methods for ordinary differential systems}, 
John Wiley, New York, 1991.

\bibitem{Lam2} J. D. Lambert; 
\emph{Computational methods in ordinary differential equations}, 
John Wiley, New York, 1973.

\bibitem{CONG} H. S. Nguyen, R. B. Sidje and N. H. Cong; 
\emph{Analysis of trigonometric implicit Runge-Kutta methods}, 
J. Comput. Appl. Math. 198 (2007) 187-207.

\bibitem{SHW} L. F. Shampine and H. A. Watts; 
\emph{Block implicit one-step methods}, Math. Comp. 23 (1969) 731-740.

\bibitem{SIMOS} T. E. Simos; 
\emph{An exponentially-fitted Runge-Kutta method for the numerical 
integration of initial-value problems with periodic or oscillating solutions},
 Comput. Phys. Commun. 115 (1998) 1-8.

\bibitem{SJ} B. P. Sommeijer; 
\emph{Explicit, high-order Runge-Kutta-Nystr\"{o}m methods for parallel computers},
 Appl. Numer. Math. 13 (1993) 221-240.

\bibitem{VAN} G. Vanden, L. Gr. Ixaru, and M. van Daele; 
\emph{Optimal implicit exponentially-fitted  Runge-Kutta}, Comput. Phys. 
Commun. 140 (2001) 346-357.

\bibitem{RAMOS} J. Vigo-Aguiar, and H. Ramos; 
\emph{Dissipative Chebyshev exponential-fitted methods for numerical 
solution of second-order differential equations}, J. Comput. Appl. Math.,
 158 (2003)  187-211.

\end{thebibliography}

\end{document}