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

\AtBeginDocument{{\noindent\small
Tenth MSU Conference on Differential Equations and Computational
Simulations. \newline
\emph{Electronic Journal of Differential Equations},
Conference 23 (2016),  pp. 1--7.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{9mm}}

\begin{document} \setcounter{page}{1}
\title[\hfilneg EJDE-2016/Conf/23 \hfil Rational logarithmic basis functions]
{Using rational logarithmic basis functions to solve
 singular differential equations}

\author[J. J. Garwood, S. N. Jator \hfil EJDE-2016/conf/23 \hfilneg]
{John J. Garwood, Samuel N. Jator}

\address{John J. Garwood \newline
Department of Mathematics and Statistics,
Austin Peay State University,
Clarksville, TN 37044, USA}
\email{jgarwood1@my.apsu.edu}

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


\thanks{Published March 21, 2016}
\subjclass[2010]{65L05, 65L06}
\keywords{Singular initial value problem; rate of convergence;
\hfill\break\indent  rational logarithmic basis function}

\begin{abstract}
 Numerical methods based on polynomial approximation perform poorly when
 applied to singular initial value problems (IVPs).  Hence, we are
 motivated to derive and implement numerical methods involving non-polynomial
 basis functions such as logarithmic and rational functions. Specifically, by
 imbedding a constant parameter into the logarithmic function, we are able to
 improve any discontinuity issues with the natural logarithm approximant.
  An efficient method is developed using the Taylor Series expansion to optimize
 the imbedded parameter.  Numerical experiments performed show that this method
 is more accurate than the improved Euler's method.  This method is implemented
 as a predictor-corrector method.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks


\section{Introduction}

We consider the initial value problem (IVP)
\begin{equation} \label{e1}
\begin{gathered}
y' = f( x,y),\\
y(a) = y_{0}, 
\end{gathered}
\end{equation}
where $y, f \in\mathbb{R}^{m}$, $x \in[a, b]$, $a, b \in\mathbb{R}$.

Most numerical methods proposed in the literature for solving IVPs 
are based on polynomial approximation. Nevertheless, methods constructed 
using polynomial basis functions perform poorly when applied to singular IVPs. 
 We define singular IVPs as those whose solution contains a point or points 
of discontinuity. Hence, several techniques involving non-polynomial basis 
functions have been proposed to cope with singular IVPs 
(see \cite{FTA,MNO,LS,Van,OMR,Ramos}).

We propose a new method constructed using rational-logarithmic basis functions 
which outperforms the Improved Euler Method (IEM).  Specifically, by imbedding 
a constant into the logarithmic function, we can improve the solutions close 
to the point of discontinuity inherent to singular IVPs. An optimized  
parameter imbedded in the method is recovered from the principal term of 
the Local Truncation Error (LTE). We then compare the performance of the method 
with the performance of the IEM.

The paper is organized as follows. In Section 2 we will examine the method's 
implementation. In Section 3 we will derive the method.  In Section 4 we will 
display numerical examples.  In Section 5 we will discuss the behavior of 
the imbedded parameter, $\alpha$, around the point of discontinuity.  
In Section 6 we will discuss opportunities for further research.

\section{Implementation of method}

The method is implemented as a predictor-corrector in a step-by-step fashion 
on the partition $H_{N}$.  The approximation at
$x_n$ is used to obtain the approximation at $x_{n+1}$, such that
\[
H _{N}:a=x_{0}<x_{1}<\dots <x_{N}=b, x_n=x_0+nh , n =0, \dots, N
\]
where $h=(b-a)/N$ is the constant step-size,  $N$ is a positive integer, 
and $n$ is the grid index.
We note that $y_n$ is the numerical approximation to the analytical
solution $y(x_n)$, and $f_n=f(x_n, y_n)$ is supplied by the differential
equation.

\section{Derivation of the method}

This method is constructed using rational-logarithmic basis functions rather 
than polynomial basis functions.

\subsection{Derivation of the corrector}
The corrector is derived by assuming that on the interval $[x_n, x_{n+1}]$,
the exact solution is approximated by the function
\begin{equation} \label{e2}
u(x) =\frac{1}{a_2(\ln[x+\alpha])^2+a_1\ln[x+\alpha]+a_0}
\end{equation}
 where $a_0, a_1, a_2$, are coefficients that are uniquely determined and 
$\alpha$ is an imbedded parameter that controls the behavior of the method. 
 We note that $\ln[x+\alpha] \in \mathbb{C}$ provided $x \neq -\alpha$.  
Thus we impose that $x\neq -\alpha$ and it follows that $u(x)$ is well defined. 
To determine the coefficients in \eqref{e2}, we require that the following 
conditions must be satisfied
\begin{gather*}
u(x_n) =  y_n,  \\
u'(x_n) =  f_n,  \\
u'(x_{n+1}) =  f_{n+1},
\end{gather*}
 This leads to a system of three equations and three unknowns, which is solved 
to obtain the coefficients in \eqref{e2}. Our method is then obtained by
 evaluating $u(x)$ at $x_{n+1}$ leading to the iterative equation
 \begin{equation} \label{e3}
y_{n+1}=\frac{y_n^2 \big(2 y_{n+1}-f_{n+1} (h+\alpha_n ) 
\ln (\frac{\alpha _n}{h+\alpha _n})\big)}
{y_{n+1} \big(f_n \alpha_n  \ln (\frac{\alpha _n}{h+\alpha _n})+2 y_n\big)}.
 \end{equation}

 Note in equation this equation, $\ln(\frac{\alpha_n}{h+\alpha_n}) \in \mathbb{R}$ 
provided $0<\alpha_n$ or $\alpha_n < -h$.  Thus by closure if for all 
$x_n, y_n,f_n,f_{n+1} \in \mathbb{R}$, then $y_{n+1} \in \mathbb{R}$.


\subsection{Derivation of the predictor}

The predictor is derived by assuming that on the interval $[x_n, x_{n+1}]$,
the exact solution is approximated by the function
\begin{equation} \label{e4}
u(x) =\frac{1}{a_2(\ln[x+\alpha])^2+a_1\ln[x+\alpha]+a_0},
\end{equation}
where $a_0, a_1, a_2$, are coefficients that are uniquely determined and 
$\alpha$ is an imbedded parameter that controls the behavior of the method. 
 Following the same assumptions for $\alpha$ as for the corrector the predictor 
is well defined. In order to determine the coefficients in \eqref{e5}, we demand
that the following conditions must be satisfied
\begin{gather*}
u(x_n) =   y_n,  \\
u'(x_n) =   f_n,  \\
u''(x_n) =  g_n,
\end{gather*}
where $g_n=\frac{df(x,y(x))}{dx}\big|_{x_n}$
This leads to a system of three equations and three unknowns, which is solved 
to obtain the coefficients in \eqref{e5}. The predictor is then obtained by evaluating
 $u(x)$ at $x_{n+1}$ to give
\begin{equation} \label{e5}
\begin{aligned}
y_{n+1}&=2 y_n^3\Big/\Big(-\alpha _n y_n \ln (\frac{\alpha _n}{h+\alpha _n}) 
\big[f_n (\ln (\frac{\alpha _n}{h+\alpha _n})-2)
 +g_n \alpha _n \ln (\frac{\alpha _n}{h+\alpha _n})\big]\\
&\quad +2 f_n^2 \alpha _n^2 \big(\ln (\frac{\alpha _n}{h+\alpha _n})\big)^2
 +2 y_n^2\Big)
\end{aligned}
\end{equation}


\subsection{Error analysis -  local truncation error}

 We carry out a Taylor series expansion about the point $x_n$ to determine 
the local truncation error as follows:
\begin{equation} \label{e6}
\begin{aligned}
LTE (y(x_n),h) 
&= y(x_n+h)-\frac{y_n^2 \big(2 y_{n+1}-f_{n+1} (h+\alpha_n ) 
\ln (\frac{\alpha _n}{h+\alpha _n})\big)}
{y_{n+1} \big(f_n \alpha_n  \ln (\frac{\alpha _n}{h+\alpha _n})+2 y_n\big)} \\
&=Ah^3 + \mathcal{O}(h^4).
\end{aligned}
\end{equation}
 This yields a principal truncation error of $Ah^3$, where
\begin{equation} \label{e7}
\begin{aligned}
A &=-\frac{ y(x_n)^2 y^{(3)}(x_n)}{6}-\frac{y(x_n)^2 y''(x_n)}{2 \alpha }
-\frac{y(x_n)^2 y'(x_n)}{6 \alpha ^2}\\
&\quad +\frac{y(x_n) y'(x_n)^2}{\alpha }
-y'(x_n)^3+y(x_n) y'(x_n) y''(x_n)
\end{aligned}
\end{equation}
Thus, the method is order $p = 2$.

\subsection{Optimization of $\alpha$}

To optimize $\alpha$, we utilize the principal truncation error.  
By setting $A = 0$ and solving for $\alpha$ we obtain
\begin{equation*} %\label{e8}
\alpha_n=\frac{-3 y_n^2 y_n''+6 y_n (y_n')^2\pm\sqrt{9 y_n^4 (y_n'')^2+12 y_n^2
 (y_n')^4-4 y_n^{(3)} y^4 y_n'-12 y_n^3 (y_n')^2 y_n''}}
{2 \big(y_n^2 y_n^{(3)}+6 (y_n')^3-6 y_n y_n'y_n''\big)},
\end{equation*}
where $\alpha_n$ is the value of $\alpha(x_n)$.  This method yields two
 possible solutions for $\alpha_n$. Taking the harmonic mean of the two terms
 yields 
\begin{equation} \label{e9}
\alpha _n=\frac{2 f_n y_n}{6 f_n^2-3 g_n y_n}.
\end{equation}
Note the arithmetic mean can be calculated as $(a+b)/2$ and the harmonic
 mean, $2ab/(a+b)$, is the reciprocal of the arithmetic mean multiplied by $ab$.  
The optimization of the Rational-Logarithmic Method (RLM) is found via the 
reciprocal of the optimization of the logarithmic method, which utilizes the 
arithmetic mean.

\section{Numerical Examples}

In this section, we give numerical examples to illustrate the accuracy of
the method.  All computations were carried out using our written code in
 Mathematica 10.0. For the sake of brevity we will not include all evaluated 
points, but will include many to show the overall trend of the approximation.

 Let $y(x_n)$ be the exact solution and $y_n$ the approximate solution on the 
partition $H_{N}$.  We find the error of the approximate solution as
 $|y(x_n)-y_n|$.  We define $E_N=|y(x_n)-y_n|$ as the maximum absolute error on 
$H_N$, and $E_n$ as the error at $x_n$. We calculate the rate of convergence 
(ROC) using the formula $\mathrm{ROC} = \log_{2}(E_{N}/E_{2N})$.


\begin{example}\label{example1} \rm
 We consider the non-singular IVP 
\begin{gather*}
y'=\frac{1}{2}-x+2y, \\
y(0)=1,
\end{gather*}
where $x\in [0,1]$ and $N=20$.
The exact solution is
\[
y(x) = e^{2x}+\frac{x}{2}
\]
 \end{example}

In our first example we will compare the performance of the rational-logarithmic 
method (RLM) with the improved Euler method (IEM).  
This is a continuous example with no discontinuous points.  
Since both methods are of order $p=2$, we may consider this a fair comparison.  
Table \ref{table1} highlights the results of this simulation.  It contains the values 
of $\alpha$, the approximate solutions given by both methods, the exact value, 
and the error at each step of the approximation.  Notice that the improved euler 
method (IEM) initially outperforms this method, but as $x \to 1$ the method becomes 
more accurate.  In this example the difference in accuracy between the two methods 
is marginal and there is no clear advantage in this method as compared with the 
improved Euler method (IEM).

 Since Example 4.1 is continuous on $[0,1]$ we may consider the rate of 
convergence of the rational-logarithmic method (RLM).  
In Table \ref{table2} we compare the rate of convergence of this method with that of 
improved Euler method (IEM).  Notice that this method converges to two confirming 
that it is of order $p=2$.  Though both converge to two it is clear from the 
table that this method converges more quickly than the improved Euler method (IEM). 
 Also this method's convergence exceeds two and then converges from the right 
instead of the left.

\begin{table}[ht]
\begin{center} {\footnotesize
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|ccccccc|}
 \hline
  $x_n $& $\alpha _n $& RLM & IEM & Exact & RLM  $E_n$  & IEM $E_n$ \\
  \hline
 0 & 0.196078 & 1 & 1 & 1 & 0 & 0 \\
 0.1 & 0.224511 & 1.27088 & 1.27103 & 1.2714 & 5.21362$\times 10^{-4}$ & 3.77758$\times 10^{-4}$ \\
 0.2 & 0.250209 & 1.59069 & 1.5909 & 1.59182 & 1.13924$\times 10^{-3}$ & 9.22647$\times 10^{-4}$ \\
 0.3 & 0.272716 & 1.97021 & 1.97043 & 1.97212 & 1.90678$\times 10^{-3}$ & 1.69012$\times 10^{-3}$ \\
 0.4 & 0.291788 & 2.42266 & 2.42279 & 2.42554 & 2.88351$\times 10^{-3}$ & 2.752$\times 10^{-3}$ \\
 0.5 & 0.307406 & 2.96414 & 2.96408 & 2.96828 & 4.14095$\times 10^{-3}$ & 4.20098$\times 10^{-3}$ \\
 0.6 & 0.319753 & 3.61435 & 3.61396 & 3.62012 & 5.76783$\times 10^{-3}$ & 6.15636$\times 10^{-3}$ \\
 0.7 & 0.329155 & 4.39732 & 4.39643 & 4.4052 & 7.87584$\times 10^{-3}$ & 8.77127$\times 10^{-3}$ \\
 0.8 & 0.336024 & 5.34243 & 5.34079 & 5.35303 & 1.06064$\times 10^{-2}$ & 1.22418$\times 10^{-2}$ \\
 0.9 & 0.340797 & 6.48551 & 6.48283 & 6.49965 & 1.41389$\times 10^{-2}$ & 1.68186$\times 10^{-2}$ \\
 1. & 0.343899 & 7.87036 & 7.86623 & 7.88906 & 1.87011$\times 10^{-2}$ & 2.28213$\times 10^{-2}$ \\
\hline
 \end{tabular} }
\end{center}
\caption{Results for example \ref{example1}}
\label{table1}
\end{table}

 \begin{example}\label{example2} \rm
We consider the singular IVP 
\begin{gather*}
y'=1+y^2, \\
y(0)=1,
\end{gather*}
where $x\in [0,1]$ and $N=20$.
The exact solution is
\[
y(x)=\tan\big(x+\frac{\pi}{4}\big)\,.
\]
 \end{example}

In Example \ref{example2} we consider a singular example.  It follows that 
$x \neq \pi/4$ since $\tan (\pi/2)$ is undefined.  For this example we will 
compare the rational-logarithmic method (RLM) with both the improved Euler method 
(IEM) and  Runge-Kutta of order 4 (RK4).  In Table \ref{table3} we see the rational-logarithmic
 Method's (RLM) performance separate from the other methods.  
In Table \ref{table4} we see only the error of each method's approximation. 
 We note that until the discontinuous point this method outperforms IEM by
 at least an order of magnitude.  RK4 outperforms this method until the 
point of discontinuity.  However, when $x_n > \frac{\pi}{2}$ this method 
is able to accurately step over the discontinuity with $E_n \approx 10^{-1}$. 
 This is a clear advantage over both of the other methods.  This is
 because the other two methods diverge at the point of discontinuity.  
This method not only outperforms another method of order $p=2$, but furthermore 
it outperforms a method of order $p=4$.

\begin{table}[htb]
\begin{center} {\footnotesize
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|ccccccc|}
 \hline
  $N$& IEM  &$E_N$  &ROC&RLM &$E_N$&ROC \\
\hline
20& 7.86623 & 2.28213$\times 10^{-2}$  &-&7.87036& 1.87011$\times 10^{-2}$&-\\
40& 7.88313&5.92887$\times 10^{-3}$& 1.94455&7.88435&  4.71061$\times 10^{-3}$&1.98914\\
80& 7.88755&1.51066$\times 10^{-3}$&1.97258&7.88788&  1.17902$\times 10^{-3}$&1.99832\\
160 &7.88867 & 3.81247$\times 10^{-4}$  &1.98638&7.88876&  2.94705$\times 10^{-4}$&2.00025\\
320&7.88896& 9.57612$\times 10^{-5}$ &1.99321&7.88898& 7.36551$\times 10^{-5}$&2.00042\\
640&7.88903& 2.39966$\times 10^{-5}$ &1.99661&7.88904& 1.84102$\times 10^{-5}$&2.00028\\
\hline
\end{tabular} }
\end{center}
\caption{Rate of convergence for example \ref{example1} }
\label{table2}
\end{table}

\begin{table}[htb]
\begin{center} {\footnotesize
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|ccccc|}
 \hline
  $x_n $& $\alpha _n $& $\text{RLM}$ &  $\text{Exact}$ & $\text{Method }E_n $ \\
 \hline
 0 & 0.333333 & 1 & 1 & 0 \\
 0.1 & 0.407675 & 1.22303 & 1.22305 & 2.36822$\times 10^{-5}$ \\
 0.2 & 0.50283 & 1.50849 & 1.5085 & 7.68901$\times 10^{-6}$ \\
 0.3 & 0.631945 & 1.89584 & 1.89577 & 7.13438$\times 10^{-5}$ \\
 0.4 & 0.821748 & 2.46525 & 2.46496 & 2.82619$\times 10^{-4}$ \\
 0.5 & 1.13636 & 3.40907 & 3.40822 & 8.50248$\times 10^{-4}$ \\
 0.6 & 1.77821 & 5.33462 & 5.33186 & 2.76644$\times 10^{-3}$ \\
 0.65 & 2.44873 & 7.3462 & 7.34044 & 5.76176$\times 10^{-3}$ \\
 0.7 & 3.89891 & 11.6967 & 11.6814 & 1.53674$\times 10^{-2}$ \\
 0.75 & 9.44124 & 28.3237 & 28.2383 & 8.54765$\times 10^{-2}$ \\
 0.8 & -22.8622 & -68.5866 & -68.4797 & 1.06918$\times 10^{-1}$ \\
 0.85 & -5.14964 & -15.4489 & -15.4579 & 8.96374$\times 10^{-3}$ \\
 0.9 & -2.89397 & -8.6819 & -8.68763 & 5.72979$\times 10^{-3}$ \\
 0.95 & -2.00545 & -6.01636 & -6.0203 & 3.9352$\times 10^{-3}$ \\
 1. & -1.52837 & -4.58511 & -4.58804 & 2.92595$\times 10^{-3}$ \\
\hline
\end{tabular} }
\end{center}
\caption{Results for example \ref{example2}}
\label{table3}
\end{table}

\begin{table}[htb]
\begin{center} {\footnotesize
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|cccc|}
 \hline
  $x_n $&  $\text{RLM }E_n$ &  $\text{IEM }E_n$ & $\text{RK4 }E_n $ \\
  \hline
 0 & 0 & 0 & 0 \\
 0.1 & 2.36822$\times 10^{-5}$  & 2.73043$\times 10^{-4}$  & 2.15331$\times 10^{-8}$  \\
 0.2 & 7.68901$\times 10^{-6}$  & 9.50794$\times 10^{-4}$  & 2.79710$\times 10^{-8}$  \\
 0.3 & 7.13438$\times 10^{-5}$  & 2.68231$\times 10^{-3}$  & 5.21260$\times 10^{-7}$ \\
 0.4 & 2.82619$\times 10^{-4}$  & 7.56687$\times 10^{-3}$  & 3.63068$\times 10^{-6}$ \\
 0.5 & 8.50248$\times 10^{-4}$  &2.42079$\times 10^{-2}$  & 2.59756$\times 10^{-5}$  \\
 0.6 & 2.76644$\times 10^{-3}$  & 1.05272$\times 10^{-1}$  & 2.91448$\times 10^{-4}$  \\
 0.65 & 5.76176$\times 10^{-3}$  & 2.77873$\times 10^{-1}$  & 1.46124 $\times 10^{-3}$ \\
 0.7 & 1.53674$\times 10^{-2}$  & 1.01467 &1.33594$\times 10^{-2}$  \\
 0.75 & 8.54765$\times 10^{-2}$   & 7.94846 & 5.4355$\times 10^{-1}$  \\
 0.8 & 1.06918$\times 10^{-1}$  &140.98 & 1392.15 \\
 0.85 & 8.96374$\times 10^{-3}$  & 3031.14  & 1.39964$\times 10^{26}$   \\
 0.9 & 5.72979$\times 10^{-3}$  &5.23821$\times 10^{9}$  & 2.69341$\times 10^{394}$ \\
 0.95 & 3.93520$\times 10^{-3}$   & 4.70556$\times 10^{34}$  & 9.52550 $\times 10^{6286}$  \\
 1. & 2.92595$\times 10^{-3}$  & 3.06426 $\times 10^{134}$  & 5.70492 $\times 10^{100567}$  \\
\hline
\end{tabular} }
\end{center}
\caption{Results for example \ref{example2} }
\label{table4}
\end{table}

\section{Behavior of $\alpha$ about the discontinuity}

The method's ability to step over the discontinuity is due to the imbedded 
parameter $\alpha$. By examining Figure \ref{fig1} it is clear that the method is 
able to accurately step over the point of discontinuity.  
If we note the optimized function 
\[
\alpha _n=\frac{2 f_n y_n}{6 f_n^2-3 g_n y_n}
\]
immediately prior to the point of discontinuity, the evaluation of 
$g_n$ is sufficient to change $\alpha_n$ from a positive value to a negative 
value. That is $3g_ny_n > 6f_n^2$.  This is confirmed in the $\alpha_n$ 
column of Table \ref{table3}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig1} 
\end{center}
\caption{Graph of exact solution with the approximate solutions of the 
Rational-Logarithmic Method}
\label{fig1}
\end{figure}


 Hence, the behavior or the imbedded constant $\alpha$ controls the behavior of the method about the point of discontinuity, permitting it to move from a relatively large positive function value to a relatively large negative function value.

\subsection*{Conclusion}
The rational-logarithmic method (RLM) is an efficient method which 
yields more accurate results than the improved Euler's method.  
It is able to accurately and efficiently overcome points of discontinuity 
in the resultant functions.  Also, the method is able to outperform methods 
of higher order in overcoming these discontinuous points.  
In future papers we will consider implementing this method in singular partial 
differential equations and systems of singular first order differential equations.

\subsection*{Acknowledgments}
 This research  was made possible through the Summer Undergraduate 
Research Fellowship of the Office of Undergraduate Research at Austin Peay 
State University.


\begin{thebibliography}{0}

\bibitem{FTA} S. O Fatunla;
\emph{Nonlinear multistep methods for initial value problems}. 
Computers \& Mathematics with Applications, 8 (1982), 231-239.

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

\bibitem{MNO} M. N. O. Ikhile;
\emph{Coefficients for Studying One-Step Rational Schemes for IVPs in ODEs: III.
Extrapolation Methods},  Computers \& Mathematics with Applications, 
47 (2004), 1463-1475.

\bibitem{LS}  J. D. Lambert;
\emph{Numerical Methods for Ordinary Differential Systems 1991}, England, John Wiley.

\bibitem{Van} F. D. Van Niekerk;
\emph{Rational one-step methods for initial value problems}, Computers Math. 
Applic. 16, (1988), 1035-1039.

\bibitem{OMR} M. R. Odekunle, N. D. Oye, S. O. Adee,  R. A. Ademiluyi;
\emph{A class of inverse Runge-Kutta schemes
for the numerical integration of singular problems}, Applied Mathematics 
and Computation, 158 (2004), 149-158.

\bibitem{Ramos}  H. Ramos, J. Vigo-Aguiar;
\emph{A new algorithm appropriate for solving singular and
singularly perturbed autonomous initial-value problems}, 
International Journal of Computer Mathematics 85, (2008) 603-611.

\end{thebibliography}

\end{document}