\documentclass{mc}
%%%%% journal  info -- DO NOT CHANGE %%%%%%%%%
\setcounter{page}{127}
\renewcommand\thisnumber{x}
\renewcommand\thisyear {2018}
\renewcommand\thismonth{xxx}
\renewcommand\thisvolume{23}
\renewcommand\datereceived{March 7, 2017}
\renewcommand\dateaccepted{September 7, 2017}
%%%%%%%%   end journal info   %%%%%%%%%%%%%
\usepackage[displaymath,mathlines,pagewise]{lineno}

%Rad zaprimljen dana: 2017-03-07

%Rad prihvacen dana: 2017-09-07

%%%{\tt amssymb, amsmath,amsthm, latexsym, amsfonts,amsbsyt, graphicx}.

%%%%% author macros %%%%%%%%%
% place your own macros HERE
%%%%% end %%%%%%%%%

%%%%%%%%%    already loaded packages -- DO NOT RELOAD  %%%%%
\usepackage{color}
% \usepackage{amsmath,amsthm}
% \usepackage{latexsym}
% \usepackage{amsfonts}
% \usepackage{amsbsyt}
% \usepackage{graphicx}
%%%%%%    end   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\I}{\mathrm{i}}
\newcommand{\E}{\mathrm{e}}
\newcommand{\D}{\mathrm{d}}
\def\RR{\mathbb{R}}
\def\ZZ{\mathbb{Z}}
\def\CC{\mathbb{C}}
\def\NN{\mathbb{N}}
\def\PP{\Cale P}

\begin{document}
%\begin{linenumbers}
%%%%% title and author(s):
% \markboth{Author(s)}{Short Title}
% \title{Title}

\markboth{M.\,Masjed-Jamei, G.\,V.\,Milovanovi\'{c} and A.\,H.\,Salehi Shayegan}{On weighted Adams-Bashforth rules}
\title[On weighted Adams-Bashforth rules]{On weighted Adams-Bashforth rules}

% single author:
% \author[AUTHOR]{AUTHOR}
% \address{address of AUTHOR}
% \email{{\tt email address of AUTHOR} (AUTHOR)}

%\author[D. Juki\'c]{Dragan Juki\'c\corrauth}
%\address{Department of Mathematics, University of Osijek, Trg Ljudevita Gaja 6, HR-31 000 Osijek, Croatia}
%\email{{\tt jukicd@mathos.hr} (D. Juki\'c)}




% multiple authors:
% Please mark \corrauth after the name of the corresponding author.
% different addresses:
\author[M.\,Masjed-Jamei, G.\,V.\,Milovanovi\'{c}, A.\,H.\,Salehi Shayegan]{Mohammad Masjed-Jamei\affil{1}, Gradimir V.\ Milovanovi\'c\affil{2}\comma\affil{3}\comma\corrauth and Amir Hossein  Salehi Shayegan\affil{1}}
\address{\affilnum{1}\ Department of Mathematics, K.\,N.\,Toosi University of Technology,
P.\,O.\,Box 16\,315--1618, Tehran, Iran\\   \affilnum{2}\ Serbian Academy of Sciences
and Arts, Kneza Mihaila 35, 11\,000 Beograd, Serbia\\
\affilnum{3}\ University of Ni\v{s}, Faculty of Sciences and Mathematics, P.\,O.\,Box 224, 18\,000 Ni\v{s}, Serbia}
%
%same address:
%\author[AUTHOR1, AUTHOR2 and AUTHOR3]{AUTHOR1, AUTHOR2\corrauth and AUTHOR3}
%\address{address of AUTHOR1, AUTHOR2 and AUTHOR3}
%
\emails{{\tt mmjamei@kntu.ac.ir}\,\,(M.\,Masjed-Jamei), {\tt gvm@mi.sanu.ac.rs} (G.\,V.\,Milovanovi\'{c}), {\tt ah.salehi@mail.kntu.ac.ir} (A.\,H.\, Salehi Shayegan)}
%




%%%%% Begin Abstract %%%%%%%%%%%
\begin{abstract}
One class of linear multistep methods for solving the Cauchy problems of the form $ y'=F(x,y) $,  $ y(x_{0})=y_{0} $, contains Adams-Bashforth rules of the form $y_{n+1}=y_{n}+h\sum_{i=0}^{k-1} B_i^{(k)} F(x_{n-i},y_{n-i})$,
where $\{ B_i^{(k)}\} _{i = 0}^{k - 1}$ are fixed  numbers. In this paper,  we propose an idea for a weighted type of Adams-Bashforth rules for solving the Cauchy problem for singular differential equations,
\[A(x)y'+B(x)y=G(x,y), \quad y(x_0)=y_0,\]
where $A$ and $B$ are two polynomials determining the well-known classical weight functions in the theory of orthogonal polynomials. Some numerical examples are also included.
\end{abstract}
%%%%% end %%%%%%%%%%%

%%%%% Keywords %%%%%%%%%%%
\keywords{weighted Adams-Bashforth rule, ordinary differential equation, linear multistep method, weight function}

%%%% AMS subject classifications %%%%
\ams{65L05, 33C45}

%%%% maketitle %%%%%
\maketitle

%%%% Start %%%%%%
\section{Introduction}

In this paper, we present an idea for constructing weighted Adams-Bashforth rules for solving Cauchy problems for singular differential equations.

There are two main approaches to increase the accuracy of a numerical method for ordinary non-singular differential equations. In the first approach (i.e., multistep methods), the accuracy is increased by considering previous information, while in the second one (i.e., multistage methods or more precisely Runge-Kutta methods), the accuracy is increased by approximating the solution at several internal points.

Multistep methods were originally proposed by Bashforth and Adams \cite{Bash} (see also \cite{3,Butcher2000,6}), where the approximate solution of the initial value problem
 \begin{equation}\label{1}
 \frac{dy}{dx}=F(x,y), \quad   y(x_0)=y_0,
\end{equation}
is considered as
\begin{equation}\label{2}
y_{n+1}=y_{n}+h\sum_{i=0}^{k-1} B_i^{(k)} F(x_{n-i},y_{n-i}).
\end{equation}


Many years later,  Moulton \cite{Moulton} (see also, \cite{Butcher2000,6}) developed a class of implicit multistep
methods, the so-called Adams-Moulton methods,
\begin{equation}\label{3}
y_{n+1}=y_{n}+h\sum_{i=-1}^{k-1} \alpha _i^{(k)} F(x_{n-i},y_{n-i}),
\end{equation}
which have some better characteristics than the previous ones.

The unknown coefficients $B_i^{(k)}$ and $\alpha _i^{(k)}$ in relations (\ref{2}) and (\ref{3}) are chosen in such a way that they have the highest possible accuracy order. These formulas are indeed special cases of the so-called linear multistep methods denoted by
\begin{equation*}
y_n=\sum_{j=1}^{k_1}\eta_jy_{n-j}+h\sum_{i=0}^{k_2}\gamma_iF(x_{n-i},y_{n-i}).
\end{equation*}

Other special cases of linear multistep methods were derived by
Nystr\"{o}m and Milne \cite{3,6}. The idea of Predictor-Corrector
methods was proposed by Milne \cite{6} in which $y_n$ is
predicted by the Adams-Bashforth methods and then corrected by the
Adams-Moulton methods.

It is not fair to talk about linear multistep methods without
mentioning the name of Germund Dahlquist. In 1956, he \cite{8} established
some basic concepts such as consistency, stability and convergence in
numerical methods and showed that if a numerical method is
consistent and stable, then it is necessarily convergent.

However, it should be noted that the above-mentioned methods are valid only for non-singular problems of type (\ref{1}). In other words, if equation (\ref{1}) is considered as an initial value problem on $ (a,b) $ in the form
\begin{eqnarray}\label{1.4}
A(x)y' = H(x,y),\quad y(a) = {y_0},
\end{eqnarray}
such that
\[A(a) = 0\quad\mbox{or}\quad A(b) = 0,\]
then it is no longer possible to use usual Adams-Bashforth methods or other numerical techniques. For this purpose, in this paper, we gave an idea for using a weighted Adams-Bashforth rule.

For constructing these weighted rules we use a similar procedure as in the case of  non-weighted formulas. Therefore, in Section 2, we give a short account of constructing the usual
Adams-Bashforth methods by using linear difference operators and the backward Newton interpolation formula. Such a procedure is applied in Section 3 for obtaining the weighted rules. By introducing the weighted local truncation error of such rules, we determine their order. Finally, in order to illustrate the efficiency of such weighted rules, we give some  numerical examples in Section 4.


\section{Computing the usual Adams-Bashforth methods}\label{sec.2}

In this section, we obtain the explicit forms of the coefficients $\{ B_i^{(k)}\} _{i = 0}^{k - 1}$ in (\ref{2}) using the backward Newton interpolation formula for $F(x,y)=F(x,y(x))$ at equidistant nodes $x_{n-\nu}=x_n-\nu h$, $\nu=0,1,\ldots,k-1$, and in the next section we apply such an approach in order to get the corresponding weighted type of Adams-Bashforth methods. Here we use standard linear difference operators $\nabla$ (the backward-difference operator), $E$ (the shifting operator), and  $1$ (the identity operator), defined by
\[\nabla f(x)=f(x)-f(x-h),\qquad Ef(x)=f(x+h)\quad\mbox{and}\quad 1f(x)=f(x).\]

Since $E^\lambda=(1-\nabla)^{-\lambda}$, we have
\begin{equation}\label{Enalam}
E^\lambda=\sum_{\nu=0}^{+\infty} (-1)^\nu\binom{-\lambda}{\nu}\nabla^\nu=	\sum_{\nu=0}^{+\infty}\frac{(\lambda)_\nu}{\nu!}\nabla^\nu,
\end{equation}
where
\[(\lambda)_\nu=\lambda(\lambda+1)\cdots(\lambda+\nu-1)=\frac{\Gamma(\lambda+\nu)}{\Gamma(\lambda)}\]
is  Pochhammer's symbol. Assuming $F_{n-\nu}\equiv F(x_{n-\nu},y(x_{n-\nu}))$ for $\nu=0,1,\ldots,k-1$ and taking the first $k$ terms of (\ref{Enalam}) for $x=x_n+\lambda h$ we get
\begin{equation}\label{BNIP}
F(x,y(x))=E^\lambda F_n=	\sum_{\nu=0}^{k-1}\frac{(\lambda)_\nu}{\nu!}\nabla^\nu F_n+r_k(F_n),
\end{equation}
where $r_k(F_n)$ denotes the corresponding error term. Using (\ref{BNIP}) we have
\begin{eqnarray}\label{IntP}
\sum_{\nu=0}^{k-1}\frac{(\lambda)_\nu}{\nu!}\nabla^\nu F_n
 &=& \sum_{\nu=0}^{k-1}\frac{(\lambda)_\nu}{\nu!}
   \sum_{i=0}^\nu (-1)^i\binom{\nu}{i}E^{-i}F_n\nonumber\\
&=& \sum_{i=0}^{k-1}\left(\frac{(-1)^i}{i!} \sum_{\nu=i}^{k-1}
   \frac{(\lambda)_\nu}{(\nu-i)!}\right)F_{n-i}=
    \sum_{i=0}^{k-1} C_i^{(k)}(\lambda)F_{n-i},
\end{eqnarray}
 where $\lambda=(x-x_n)/h$ and
 \begin{eqnarray}\label{Ciklam}
C_i^{(k)}(\lambda ) = \frac{{{{( - 1)}^i}}}{{i!}}\sum\limits_{\nu  = i}^{k - 1} {\frac{{{{(\lambda )}_\nu }}}{{(\nu  - i)!}}}  = {( - 1)^i}\frac{{{{(\lambda )}_i}}}{{i!}}\left( \begin{gathered}
  \lambda  + k - 1 \\
  k - 1 - i \\
\end{gathered}  \right),
\end{eqnarray}
because, based on induction, we have
\begin{equation*}
\begin{gathered}
  \sum\limits_{\nu  = i}^k {\frac{{{{\left( \lambda  \right)}_\nu }}}{{\left( {\nu  - i} \right)!}}}  = \sum\limits_{\nu  = i}^{k - 1} {\frac{{{{\left( \lambda  \right)}_\nu }}}{{\left( {\nu  - i} \right)!}}}  + \frac{{{{\left( \lambda  \right)}_k}}}{{\left( {k - i} \right)!}} \hfill \\
  \qquad\qquad\quad\, = {\left( \lambda  \right)_i}\left( {\begin{array}{*{20}{c}}
  {\lambda  + k - 1} \\
  {k - 1 - i}
\end{array}} \right) + {\left( \lambda  \right)_i}\left( {\begin{array}{*{20}{c}}
  {\lambda  + k - 1} \\
  {k - i}
\end{array}} \right) = {\left( \lambda  \right)_i}\left( {\begin{array}{*{20}{c}}
  {\lambda  + k} \\
  {k - i}
\end{array}} \right). \hfill \\
\end{gathered}
\end{equation*}

Now, integrating (\ref{1}) over $(x_n,x_{n+1})$ and approximating $F(x,y)$ by its backward Newton interpolation polynomial (\ref{IntP}) yield
\begin{equation}\label{2.5}
  y({x_{n + 1}}) - y({x_n})= \int_{{x_n}}^{{x_{n + 1}}} F (x,y)\text{d}x \approx h\sum_{i=0}^{k-1} B_i^{(k)}F_{n-i},
\end{equation}
where
\begin{equation}\label{Bik}
B_i^{(k)}=\int_0^1 C_i^{(k)}(\lambda)\text{d}\lambda=\frac{(-1)^i}{i!}\int_0^1(\lambda)_i \binom{\lambda+k-1}{k-1-i}\text{d}\lambda.
\end{equation}

Table \ref{tb1} shows the values $B_i^{(k)}$ in (\ref{2}) for $k=2,3,\ldots,8$.\\
\begin{table}[!htb]
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
$k$ & $2$ & $3$ & $4$ & $5$ & $6$ & $7$ & $8$ \\
\hline
~$B_0^{(k)}$&$\frac{3}{2}$&$\frac{23}{12}$&$\frac{55}{24}$&$\frac{1901}{720}$&$\frac{4277}{1440}$&$\frac{198721}{60480}$&$\frac{16083}{4480}$\\[2mm]
~$B_1^{(k)}$&$-\frac{1}{2}$&$-\frac{4}{3}$&$-\frac{59}{24}$&$-\frac{1387}{360}$&$-\frac{2641}{480}$&$-\frac{18637}{2520}$&$-\frac{1152169}{120960}$\\[2mm]
~$B_2^{(k)}$& &$\frac{5}{12}$&$\frac{37}{24}$&$\frac{109}{30}$&$\frac{4991}{720}$&$\frac{235183}{20160}$&$\frac{242653}{13440}$\\[2mm]
~$B_3^{(k)}$& & &$-\frac{3}{8}$&$-\frac{637}{360}$&$-\frac{3649}{720}$&$-\frac{10754}{945}$&$-\frac{296053}{13440}$\\[2mm]
~$B_4^{(k)}$& &  & &$\frac{251}{720}$&$\frac{959}{480}$&$\frac{135713}{20160}$&$\frac{2102243}{120960}$\\[2mm]
~$B_5^{(k)}$& & & & &$-\frac{95}{288}$&$-\frac{5603}{2520}$&$-\frac{115747}{13440}$\\[2mm]
~$B_6^{(k)}$& & & & & &$\frac{19087}{60480}$&$\frac{32863}{13440}$\\[2mm]
~$B_7^{(k)}$& ~&~ ~&~ ~&~ ~&~ ~&~ ~&~$-\frac{5257}{17280}$\\
\hline
\end{tabular}
\medskip\caption[]{\rm\em The coefficients of usual Adams-Bashforth formulae }\label{tb1} \end{table}


By assuming that all previous values $y_{n-i}$, $i=0,1,\ldots,k-1$, are exact, i.e., $y_{n-i}=y(x_{n-i})$, $i=0,1,\ldots,k-1$, (\ref{2.5}) gives the
$k$-step method  (\ref{2}). This $k$-step method,  known also as
 {\it the $k$th-order Adams-Bashforth method},  can be  written in the form
\[y_{n+k}-y_{n+k-1}=h\sum_{j=0}^{k-1} \beta_j^{(k)} F_{n+j},\]
where $\beta_j^{(k)}=B_{k-1-j}^{(k)}$, $j=0,1,\ldots,k-1$.

According to (\ref{BNIP}), {\it the local truncation error} of this method
at the point $x_{n+k}\in[a,b]$ can be expressed in the form
\begin{equation}\label{LTE}
(T_h)_{n+k}=\frac{y(x_{n+k})-y(x_{n+k-1})}{h}-\sum_{j=0}^{k-1} \beta_j^{(k)}y'(x_{n+j})=\int_0^1 r_k(F_{n+k-1}){\D}\lambda, 	
\end{equation}
where $x\mapsto y(x)$ is the exact solution of the Cauchy problem (\ref{1}).
If $y\in C^{k+2}[a,b]$, then  (\ref{LTE})
can be expressed as (cf. \cite[pp.~409--410]{WG})
\begin{equation}\label{Th-st}
(T_h)_{n+k}=C_k y^{(k+1)}(\xi_k)h^k=C_k y^{(k+1)}(x_n)h^k +O(h^{k+1}), 	
\end{equation}
where $x_n<\xi_k<x_{n+k-1}$.
 In the simplest case ($k=1$), we have the well-known Euler method, $y_{n+1}-y_n=hF_n$. The so-called {\it error constants} $C_k$ in {\it the main term of the local truncation error} for $k=1,2,3,4$ and $5$ are
\begin{equation}\label{errconst}
C_1=\frac{1}{2},\quad  C_2= \frac{5}{12},\quad C_3=\frac{3}{8},\quad C_4=\frac{251}{720},\quad C_5=\frac{95}{288},	
\end{equation}
respectively. Details on multistep methods, including convergence, stability and estimation of global errors  $e_n=y_n-y(x_n)$, can be found in \cite{6,WG,gvm88}.


\begin{remark}
{\rm These coefficients $B_i^{(k)}$ can  also be expressed in terms of the first kind Stirling numbers $\mathbf{S}(n,k)$, which are defined by
\[\prod\limits_{i = 0}^{n - 1} {(x - i)}  = \sum\limits_{k = 0}^n {\mathbf{S}} (n,k){x^k},\]
(see \cite{7,20,22,18}).
Namely, for each $k\in \mathbb{N}$,  coefficients \eqref{Bik} can be explicitly represented in terms of the first kind Stirling numbers as
 \begin{equation}\label{StirExp}
B_i^{(k)}=\sum_{\nu=i}^{k-1}\left(\frac{\displaystyle\sum_{j=0}^\nu\frac{(-1)^j}{j+1}\mathbf{S}(\nu,j)}{(-1)^\nu \nu!+\displaystyle\sum_{j=1}^\nu i^j(j+1)\mathbf{S}(\nu+1,j+1)}\right),\quad   i=0,1,\ldots,k-1.
\end{equation}}
\end{remark}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Weighted type of Adams-Bashforth methods}\label{sec3}

In this section, we study the Cauchy problem for a special type of differential equations of the first order given on a finite interval,  on a half line or on the real line, which can be considered, without loss of generality, as $(-1,1)$, $(0,+\infty)$, and $(-\infty,+\infty)$. Thus, we consider the following initial value problem on $(a,b)$
\begin{equation}\label{WCP}
A(x)y'+B(x)y=G(x,y), \quad y(x_0)=y_0,
\end{equation}
where $A$ and $B$ are two polynomials determining the well-known classical weight functions in the theory of orthogonal polynomials (cf. \cite[p.~122]{Mas_Mil2008}). Such polynomials and weight functions are given in Table \ref{tab2}, where $\alpha,\beta,\gamma>-1$.

\begin{table}[!htb]
\centering
\begin{tabular}{|c|c|c|c|}
\hline
$(a,b)$ & $w(x)$ & $A(x)$ & $B(x)$ \\
\hline
$(-1,1)$ & $(1-x)^\alpha(1+x)^\beta$ & $1-x^2$ & $\beta-\alpha-(\alpha+\beta+2)x$ \\ %[2mm]
$(0,+\infty)$ & $x^\gamma {\E}^{-x}$ & $x$ & $\gamma+1-x$\\%[2mm]
$(-\infty,+\infty)$ &  ${\E}^{-x^2}$ & $1$ & $-2x$\\
\hline
\end{tabular}
\medskip\caption{\rm \em Classical weight functions and corresponding polynomials $A$  and $B$}\label{tab2}\end{table}

Let again $\{x_k\}$ be a system of equidistant nodes with the step $h$, i.e., $x_k=x_0+kh \in [a,b]$.

Since the differential equation of the weight function is as (cf. \cite[p.~122]{Mas_Mil2008})
\[(Aw)'=Bw, \]
after multiplying by $w(x)$, our initial value problem (\ref{WCP}) becomes
\[A(x)w(x)y'+B(x)w(x)y=w(x)G(x,y),\quad y(x_0)=y_0,\]
which is equivalent to
\begin{equation}\label{SLWCP}
(A(x)w(x)y)'=w(x)G(x,y),\quad y(x_0)=y_0.	
\end{equation}

Now, integrating from both sides of (\ref{SLWCP}) over $ [x_{n},x_{n+1}] $ yields
\begin{equation}\label{IntegSL}
A(x_{n+1})w(x_{n+1})y(x_{n+1})-A(x_n)w(x_n)y(x_n)=\int_{x_n}^{x_{n+1}}w(x)G(x,y)\text{d}x.
\end{equation}

Let $x=x_n+\lambda h$. Similar to relations \eqref{BNIP}, \eqref{IntP} and \eqref{2.5}, the right-hand side of  \eqref{IntegSL} can be
written in the form
\begin{equation}\label{EQNIP}
\int_{x_n}^{x_{n+1}}w(x)G(x,y)\text{d}x= h\int_0^1 w(x_n+\lambda h)\left\{\sum_{\nu=0}^{k-1}\frac{(\lambda)_\nu}{\nu!}\nabla^\nu G_n	+r_k(G_n)\right\}\text{d}\lambda, 	
\end{equation}
and approximated as
\begin{eqnarray*}
\int_{x_n}^{x_{n+1}}w(x)G(x,y)\text{d}x
&\approx &h\int_0^1 w(x_n+\lambda h)\sum_{\nu=0}^{k-1}\frac{(\lambda)_\nu}{\nu!}\left(\sum_{i=0}^\nu (-1)^i\binom{\nu}{i} G_{n-i}\right)\text{d}\lambda \\
&=&h\sum_{i=0}^{k-1}\left(\frac{(-1)^i}{i!}\sum_{\nu=i}^{k-1}\frac{1}{(\nu-i)!}\int_0^1 w(x_n+\lambda h)(\lambda)_\nu \text{d}\lambda\right)G_{n-i}\\
&=&h\sum_{i=0}^{k-1} B_i^{(k)}(h,x_n) G_{n-i},	
\end{eqnarray*}
 where
 \begin{equation}\label{Bikhxn}
 B_i^{(k)}(h,x_n)=\int_0^1 w(x_n+\lambda h)C_i^{(k)}(\lambda)\text{d}\lambda,
\end{equation}
$G_{n-i}\equiv G(x_{n-i},y(x_{n-i}))$, $C_i^{(k)}(\lambda)$ is given by \eqref{Ciklam} and
$r_k(G_n)$ is  the error term in the corresponding
  backward Newton interpolation formula for $G(x,y)=G(x,y(x))$ at equidistant nodes $x_{n-i}=x_n-i h$, $i=0,1,\ldots,k-1$. Hence, the approximate form of (\ref{IntegSL}) becomes
\begin{equation}\label{IntegSLFin}
A(x_{n+1})w(x_{n+1})y(x_{n+1})-A(x_n)w(x_n)y(x_n)=h\sum_{i=0}^{k-1} B_i^{(k)}(h,x_n) G_{n-i},
\end{equation}
where the coefficients $B_i^{(k)}(h,x_n)$ depend on $h$ and $x_n$. As in the case of the standard Adams-Bashforth methods, by assuming that all previous values $y_{n-i}$, $i=0,1,\ldots,k-1$, are exact, i.e., $y_{n-i}=y(x_{n-i})$, $i=0,1,\ldots,k-1$, (\ref{IntegSLFin}) gives our {\it weighted $k$-step method}
\begin{equation}\label{IntegSLFinNM}
A(x_{n+1})w(x_{n+1})y_{n+1}-A(x_n)w(x_n)y_n=h\sum_{i=0}^{k-1} B_i^{(k)}(h,x_n) G_{n-i},\quad n\ge k-1,
\end{equation}
where $G_{n-i}\equiv G(x_{n-i},y_{n-i})$, $i=0,1,\ldots,k-1$.

The mentioned dependence of the coefficients $B_i^{(k)}(h,x_n)$  on the stepsize $h$ and $x_n$ makes these methods fundamentally different from the standard ones.

Similarly to (\ref{LTE}),  we can here define  the corresponding {\it weighted local truncation error} at the point $x_{n+k}\in[a,b]$ as
\begin{eqnarray*}
(T_h^w)_{n+k}&=&\Biggl\{\frac{1}{h}\left[A(x_{n+k})w(x_{n+k}	)y(x_{n+k})-A(x_{n+k-1})w(x_{n+k-1})y(x_{n+k-1})\right]\\
&&\qquad -\sum_{j=0}^{k-1}B_{k-j-1}^{(k)}(h,x_{n+k-1})G(x_{n+j},y(x_{n+j}))\Biggr\}\frac{1}{A(x_{n+k})w(x_{n+k}	)},
\end{eqnarray*}
where $x\mapsto y(x)$ is the exact solution of the Cauchy problem (\ref{WCP}).


Then, according to (\ref{IntegSL}), with $n:=n+k-1$,
and using (\ref{EQNIP}) and (\ref{IntegSLFin}), we obtain
\[(T_h^w)_{n+k}=\frac{1}{A(x_{n+k})w(x_{n+k})h}\int_{x_{n+k-1}}^{x_{n+k}}w(x)r_k(G_{n+k-1}){\D}x.\]
The first term omitted  in the summation on the right-hand side  in
(\ref{EQNIP}) is a good approximation  of the truncation error. We will call this  quantity  {\it  the main term of the truncation error} and denote by $(\widehat T_h^w)_{n+k}$.

\begin{proposition}\label{Prop1} Let the exact solution of the singular Cauchy problem $(\ref{WCP})$ be sufficiently smooth, as well as the function $x\mapsto g(x)=G(x,y(x))$. Then, the
 main term of the truncation error at the point
$x_{n+k}$ can be expressed in the form
\begin{equation}\label{divdf}
(\widehat T_h^w)_{n+k}=\frac{h^k g^{(k)}(\xi_k)}{A(x_{n+k})w(x_{n+k})}\int_0^1\binom{\lambda+k-1}{k}w(x_{n+k-1}+\lambda h){\D}\lambda,
\end{equation}
where
\begin{equation}\label{Izrgk}
 g^{(k)}(x)=A(x)y^{(k+1)}+\left[B(x)+kA'(x)\right]y^{(k)}+
k\Bigl[B'(x)+\frac{1}{2}(k-1)A''(x)\Bigr]y^{(k-1)}
\end{equation}
and $x_{n-1}<\xi_k<x_{n+k-1}$.
\end{proposition}

\begin{proof}
According to (\ref{EQNIP}),  we have
\[(\widehat T_h^w)_{n+k}=\frac{\nabla^k G_{n+k-1}}{A(x_{n+k})w(x_{n+k})k!}\int_0^1 w(x_{n+k-1}+\lambda h)(\lambda)_k{\D}\lambda,\]
where the factor in front of the integral can be expressed in terms of
divided differences as (cf. \cite[p.~410]{WG})
\[
\frac{\nabla^k g_{n+k-1}}{k!}=h^k\left[x_{n+k-1},x_{n+k-2},\ldots,x_{n-1}\right]g.	
\]

On the other hand, supposing that the function $x\mapsto g(x)=G(x,y(x))$ is sufficiently smooth, we can write
\[\left[x_{n+k-1},x_{n+k-2},\ldots,x_{n-1}\right]g=\frac{g^{(k)}(\xi_k)}{k!},\]
 where $\xi_k$ is between the smallest and the largest of these points.
In order to calculate these derivatives,
\[g'(x)=\frac{\partial G}{\partial x}+	\frac{\partial G}{\partial y}y',\quad
g''(x)=\frac{\partial^2 G}{\partial x^2}+\left[2\frac{\partial^2 G}{\partial x\partial y}+\frac{\partial^2 G}{\partial y^2}y'\right]y'+
\frac{\partial G}{\partial y}y'',\ \ \mbox{etc.}\]
we use the relation $g(x)=G(x,y(x))=A(x)y'(x)+B(x)y(x)$. Since $A(x)$ is a polynomial of degree at most two and $B(x)$ is a polynomial of first degree,  these derivatives  can be calculated much simpler
for each $k\ge0$ in the  form (\ref{Izrgk}).

In this way, we obtain (\ref{divdf}).
\end{proof}


Formally, (\ref{divdf}) is of the same form as
(\ref{Th-st}), i.e.,
$C_k g^{(k)}(\xi_k)h^k$, where
\[C_k=C_k(h,x_n)=\frac{1}{A(x_{n+k})w(x_{n+k})}\int_0^1\binom{\lambda+k-1}{k}w(x_{n+k-1}+\lambda h){\D}\lambda\]
and $g^{(k)}$ is given by (\ref{Izrgk}). In the case of  standard Adams-Bashforth methods, $C_k$ are the  error constants  given by (\ref{errconst}) and they are independent of the stepsize $h$ and $x_n$. Also,   instead of $g^{(k)}$, there is only the derivative $y^{(k+1)}$ in (\ref{Th-st}). Because of these  differences, the actual order of the weighted methods can be reduced (see examples in Section \ref{sec:ne}).

In the sequel, we consider the three cases given previously in Table~\ref{tab2}.

\subsection{Case $(a,b)=(-1,1)$}\label{case1}
 Consider the Cauchy problem of Jacobi type
\[(1 - {x^2})y' + \left( {\beta  - \alpha  - (\alpha  + \beta  + 2)x} \right)y = G(x,y),{\text{ }}\quad y({x_0}) = {y_0},\]
where $x_n=x_0+nh\in(-1,1)$.

In this case, relation \eqref{IntegSLFin} is reduced to
\begin{equation}\label{RecJT}
y_{n+1}= d(h,x_n)y_n+h\sum_{i=0}^{k-1} D_i^{(k)}(h,x_n) G_{n-i},	
\end{equation}
where
\[d(h,{x_n}) = \frac{{A({x_n})w({x_n})}}{{A({x_{n + 1}})w({x_{n + 1}})}} = \frac{(1-x_n)^{\alpha+1}(1+x_n)^{\beta+1}}
{(1-x_{n+1})^{\alpha+1}(1+x_{n+1})^{\beta+1}}\]
and
\[D_i^{(k)}(h,{x_n}) = \frac{1}{A(x_{n + 1})w(x_{n + 1})}\int_0^1 w(x_n+\lambda h)C_i^{(k)}(\lambda) d\lambda.\]
Putting
\[c(h,{x_n}) = \frac{(1-x_n)^{\alpha}(1+x_n)^{\beta}}
{(1-x_{n+1})^{\alpha+1}(1+x_{n+1})^{\beta+1}}\]
and
\begin{equation}\label{KoefJT}
\Phi _i^{(k)}(h,{x_n}) = \int_0^1 {{{\left( {1 - {x_n} - \lambda h} \right)}^\alpha }} {\left( {1 + {x_n} + \lambda h} \right)^\beta }C_i^{(k)}(\lambda )\text{d}\lambda,
\end{equation}
 rule \eqref{RecJT} can be simplified as
\begin{equation}\label{RecJTFin1}
y_{n+1}=c(h,x_n)\left((1-x_n^2)y_n+h\sum_{i = 0}^{k - 1}  \Phi_i^{(k)}(h, x_n) G_{n - i}\right).
\end{equation}

In the special (Legendre) case $\alpha=\beta=0$, \eqref{KoefJT} is reduced to
\begin{equation}\label{KOEF}
\Phi_i^{(k)}(h,x_n)=
\int_0^1 C_i^{(k)}(\lambda)d\lambda=B_i^{(k)},
\end{equation}
where $B_i^{(k)}$ are the same coefficients as for standard (non-weighted) Adams-Bashforth formulas given by \eqref{Bik}.

\subsection{Case $(a,b)=(0,\infty)$}\label{case2}

Now, consider the Cauchy problem of Laguerre type
\[xy'+(\gamma+1-x)y=G(x,y), \quad y(x_0)=y_0,\]
in which $x_n=nh$ for $n=0,1,\ldots$, and the main relation \eqref{IntegSLFin} is reduced to the corresponding equation \eqref{RecJT} with
\[d(h,x_n) = \frac{{A({x_n})w({x_n})}}{{A({x_{n + 1}})w({x_{n + 1}})}} = {\E}^h \left(\frac{x_n}{x_{n + 1}} \right)^{\gamma  + 1},\]
and
\[D_i^{(k)}(h,{x_n}) = \frac{{B_i^{(k)}(h,{x_n})}}{{A({x_{n + 1}})w({x_{n + 1}})}} = \frac{{{{\E}^h}}}{{x_{n + 1}^{\gamma  + 1}}}\int_0^1 {{{({x_n} + \lambda h)}^\gamma }} {{\E}^{ - \lambda h}}C_i^{(k)}(\lambda )\text{d}\lambda .\]
In other words, we have
\[D_i^{(k)}(h,{x_n}) = x_n^{ - (\gamma  + 1)}d(h,{x_n})\Phi _i^{(k)}(h,{x_n}),\]
such that
\begin{equation}\label{Fix_n}
\Phi _i^{(k)}(h,{x_n}) = \int_0^1 {{{\left( {{x_n} + \lambda h} \right)}^\gamma }{{\E}^{ - \lambda h}}C_i^{(k)}(\lambda )\text{d}\lambda .}
\end{equation}
Hence, the relation corresponding to \eqref{IntegSLFin} takes the form
\begin{equation}\label{RecJTFin}
{y_{n + 1}} = {{\E}^h}{\left( {\frac{{{x_n}}}{{{x_{n + 1}}}}} \right)^{\gamma  + 1}}{y_n} + \frac{{{{\E}^h}h}}{{x_{n + 1}^{\gamma  + 1}}}\sum\limits_{i = 0}^{k - 1} {\Phi _i^{(k)}} (h,{x_n}){G_{n - i}}.
\end{equation}
Note that when $\gamma=0$, the coefficients \eqref{Fix_n} are independent of  $x_n$ and
\begin{equation}\label{KOEFExp}
\Phi _i^{(k)}(h,{x_n}) = \Phi _i^{(k)}(h) = \int_0^1 {{{\E}^{ - \lambda h}}C_i^{(k)}(\lambda ){\text{d}}\lambda .}
\end{equation}
For instance, for $k=1(1)5$ and $i=0,1,\ldots,k-1$, relation (\ref{KOEFExp})   gives
\begin{eqnarray*}
&&\underline{k=1}:\\
&&\Phi_0^{(1)} (h)=\frac{1-{\E}^{-h}}{h};\\[2mm]
&&\underline{k=2}:\\
&&\Phi_0^{(2)} (h) =\frac{1+h-{\E}^{-h} (1+2 h)}{h^2},\\
&&\Phi_1^{(2)} (h) =-\frac{1-{\E}^{-h}(1+h)}{h^2};\\[2mm]
&&\underline{k=3}:\\
&&\Phi_0^{(3)}(h) = \frac{2 + 3 h + 2 h^2-{\E}^{-h}(2+5h+6h^2)}{2h^3},\\
&&\Phi_1^{(3)}(h) = -\frac{2(1+h)-{\E}^{-h}(2 + 4 h + 3 h^2)}{h^3},
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\
&&\Phi_2^{(3)}(h) = \frac{2 + h -{\E}^{-h}(2+3h+2h^2)}{2h^3};\\[2mm]
% \end{eqnarray*}
% \begin{eqnarray*}
&&\underline{k=4}:\\
&& \Phi_0^{(4)} (h) =\frac{ 6+12 h+11 h^2+6 h^3-2{\E}^{-h}(3+9h+13 h^2+12 h^3)}{6 h^4},\\	
&&\Phi_1^{(4)} (h) =-\frac{2(3+5h+3h^2)-{\E}^{-h}(6+16h+19h^2+12h^3)}{2h^4},\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\
&&\Phi_2^{(4)} (h) =\frac{6 + 8 h + 3 h^2-2{\E}^{-h}(3 + 7 h + 7 h^2 + 4 h^3)}{2h^4},\\
&&\Phi_3^{(4)} (h) =-\frac{2(3 + 3 h + h^2)- {\E}^{-h}(6 + 12 h + 11 h^2 + 6 h^3)}{6 h^4}; \\[2mm]
%\end{eqnarray*}
%\begin{eqnarray*}
&&\underline{k=5}:\\
&& \Phi_0^{(5)} (h) =\frac{12 + 30 h + 35 h^2 + 25 h^3 + 12 h^4-{\E}^{-h}(12 + 42 h + 71 h^2 + 77 h^3 + 60 h^4)}{12 h^5},\\	
&&\Phi_1^{(5)} (h) =-\frac{2(12 + 27 h + 26 h^2 + 12 h^3)-{\E}^{-h}(24 + 78 h + 118 h^2 + 107 h^3 + 60 h^4)}{6h^5},\\
&&\Phi_2^{(5)} (h) =\frac{12 + 24 h + 19 h^2 + 6 h^3-{\E}^{-h}(12 + 36 h + 49 h^2+ 39 h^3 + 20 h^4)}{2h^5},\\
&&\Phi_3^{(5)} (h) =-\frac{2(12 + 21 h + 14 h^2 + 4 h^3)- {\E}^{-h}(24 + 66 h + 82 h^2 + 61 h^3 + 30 h^4)}{6 h^5},\\
&&\Phi_4^{(5)} (h) =\frac{12 + 18 h + 11 h^2 + 3 h^3- {\E}^{-h}(12 + 30 h + 35 h^2 + 25 h^3 + 12 h^4)}{12 h^5}.
\end{eqnarray*}

\begin{remark} {\rm According to (\ref{KOEFExp}), it is clear that
$\lim\limits_{h\to 0}\Phi_i^{(k)} (h)=B_i^{(k)}$, $i=0,1,\ldots,k-1$,
where $B_i^{(k)}$ are the  coefficients of standard (non-weighted) Adams-Bashforth formulas given by \eqref{Bik}.
}	
\end{remark}


\begin{remark} {\rm Consider the Cauchy-Laguerre problem $xy'+(1-x)y=y$  (i.e., $G(x,y)=y)$,
which is simplified as
\[y'=y,\quad y(0)=1,\]
with the exact solution $y={\E}^x$. Considering the simplest method (for  $k=1$) gives
\[{y_{n + 1}} = \frac{{{{\text{e}}^h}}}{{1 + \frac{h}{{{x_n}}}}}\left( {{y_n} + \frac{h}{{{x_n}}}\Phi _0^{(1)}(h){G_n}} \right) = \frac{{{{\text{e}}^h}}}{{{x_{n + 1}}}}\left( {{x_n}{y_n} + (1 - {{\text{e}}^{ - h}})G({x_n},{y_n})} \right).\]


Now, substituting $x_n=nh$ in the above relation yields
\[{y_{n + 1}} = \frac{(1 + nh){\E}^h - 1}{(n + 1)h}y_n,\quad \mbox{with}\quad {y_0} = 1.\]
For example, we have
\begin{eqnarray*}
&&y_1=\frac{{\E}^h-1}{h},\quad  y_2=\frac{(1+h){\E}^h-1}{2h}\cdot \frac{{\E}^h-1}{h},\\[2mm]
&&y_3=\frac{(1+2h){\E}^h-1}{3h}\cdot \frac{(1+h){\E}^h-1}{2h}\cdot \frac{{\E}^h-1}{h},	
\end{eqnarray*}
or, in general,
\[ y_n = \prod_{\nu  = 1}^n\frac{[1 + (\nu -1)h]{\E}^h-1}{\nu h}.\]

The method is convergent, i.e.,
\begin{equation}\label{konver}
\lim_{\scriptsize\begin{array}{c}
n\to+\infty\\
(nh=x=\mbox{const})	
\end{array}} y_n=
\lim_{n\to+\infty}\prod_{\nu=1}^n \frac{\left[1+(\nu-1)\dfrac{x}{n}\right]{\E}^{x/n}-1}{\dfrac{\nu x}{n}}={\E}^x.	
\end{equation}

In order to prove (\ref{konver}) we define a sequence
$\{a_n\}_{n\in\NN}$ by
\[a_n=\sum_{\nu=1}^n \log\left\{\frac{\left[1+(\nu-1)\dfrac{x}{n}\right]{\E}^{x/n}-1}{\dfrac{\nu x}{n}}\right\}^n\]
and apply the well known Stolz-Ces\`aro theorem. Namely, if we prove the convergence
\[\lim_{n\to+\infty}\frac{a_n-a_{n-1}}{n -(n-1)}=\lim_{n\to+\infty}
\log \left\{\frac{\left[1+\left(1-\dfrac{1}{n}\right)x\right]{\E}^{x/n}-1}{x}\right\}^n=L,\]
then  the limit $\lim\limits_{n\to+\infty}\dfrac{a_n}{n}$ also exists and it is equal to $L$.
Since
\[\frac{1}{x}\left\{\left[1+\left(1-\dfrac{1}{n}\right)x\right]{\E}^{x/n}-1\right\}=1+\frac{x}{n}+O(n^{-2}), \]
we conclude that $L=\log{\E}^x=x$. Therefore, we obtain (\ref{konver}), because
\begin{eqnarray*}
\lim_{n\to+\infty}\frac{a_n}{n}&=&\lim_{n\to+\infty}\sum_{\nu=1}^n
\log\left\{\frac{\left[1+(\nu-1)\dfrac{x}{n}\right]{\E}^{x/n}-1}{\dfrac{\nu x}{n}}\right\}\\
&=&\lim_{n\to+\infty}\log\left\{\prod_{\nu=1}^n \frac{\left[1+(\nu-1)\dfrac{x}{n}\right]{\E}^{x/n}-1}{\dfrac{\nu x}{n}}\right\}=x.	
\end{eqnarray*}}
\end{remark}


\subsection{Case $(a,b)=(-\infty,\infty)$} \label{case3}

Now, consider the Cauchy problem of Hermite type
\[y'-2xy=G(x,y),\quad y(x_0)=y_0,\]
in which $x_n=x_0+nh\in(-\infty,\infty)$ and the main relation \eqref{IntegSLFin} is reduced to the corresponding
equation \eqref{RecJT} with
\[
d(h,x_n)=\frac{{\E}^{-x_n^2}}{{\E}^{-x_{n+1}^2}}=
{\E}^{x_{n+1}^2-x_n^2}={\E}^{h(2x_n+h)},	
\]
and
\[D_i^{(k)}(h,{x_n}) = \frac{{B_i^{(k)}(h,{x_n})}}{{{{\text{e}}^{ - x_{n + 1}^2}}}} = {{\text{e}}^{x_{n + 1}^2}}\int_0^1 {{{\text{e}}^{ - {{({x_n} + \lambda h)}^2}}}} C_i^{(k)}(\lambda ){\text{d}}\lambda .\]
In other words, we have
$D_i^{(k)}(h,{x_n}) = d(h,{x_n})\Phi _i^{(k)}(h,{x_n})$,
such that
\[\Phi _i^{(k)}(h,{x_n}) = \int_0^1 {{{\text{e}}^{ - (2{x_n}\lambda h + {\lambda ^2}{h^2})}}} C_i^{(k)}(\lambda ){\text{d}}\lambda.\]
Hence, the relation corresponding to \eqref{IntegSLFin} takes the form
\begin{equation}\label{Hermite}
{y_{n + 1}} = {{\text{e}}^{h(2{x_n} + h)}}\left( {{y_n} + \sum\limits_{i = 0}^{k - 1} {\Phi _i^{(k)}} (h,{x_n}){G_{n - i}}} \right).
\end{equation}

\begin{remark}\label{stvr}
 {\rm As in the case of non-weighted methods, in applications of these methods for $k>1$, we need the additional starting values $y_i=s_i(h)$, $i=1,\ldots,k-1$, such that $\lim\limits_{h\to 0}=y_0$ (cf. \cite[pp.~32--36]{gvm88}).}	
\end{remark}


\section{Numerical examples}\label{sec:ne}

In order to illustrate the efficiency of our method, in this section we give two numerical examples for singular Cauchy problems on $(0,\infty)$ and $(-1,1)$. In particular, the weighted Adams-Bashforth methods with respect to the standard Laguerre weight given in Subsection 3.2 have the simplest form and they can find adequate application in solving weighted singular Cauchy problems. The third case when $(a,b)=(-\infty,\infty)$ is not interesting for applications because  equation (\ref{WCP}) is not singular in a finite domain.

\begin{example}\label{ex1} {\rm We first consider a singular Cauchy problem
\[(1-x^2)y'-3x y=\frac{y^2 \left(\left(1-x^2\right) \tan (x)+4 x+1\right) \sec (x)}{x-1},\quad y(-1)=2\cos1.\]
Here,
\[G(x,y)=\frac{y^2 \left(\left(1-x^2\right) \tan (x)+4 x+1\right) \sec (x)}{x-1}+x y\]
and the exact solution of this  problem is given by
$y(x)=(1-x)\cos x$.

In order to solve the problem for $x\in[-1,0]$, we apply the $k$-step method (\ref{RecJTFin1}), with  $\alpha=\beta=0$ (Legendre case) and $\Phi_i^{(k)}(h,x_n)$ given by (\ref{KOEF}). For the sake of simplicity, in the case  $k>1$,  for starting values we use the exact values. Otherwise, some other ways must be applied (see Remark \ref{stvr}).

Relative errors obtained by this $k$-step method when $k=1,2,\ldots,5$,  for $h=0.02$ and $h=0.01$, are displayed in $\log$-scale in Figures \ref{fig1} and \ref{fig1a}, respectively.
%\medskip
\begin{figure}[h]
\centering
\includegraphics[width=10.cm]{fig1.eps}
%\includegraphics[width=10.cm]{figure1.eps}
\medskip
\caption{\rm\em Relative errors for methods (\ref{RecJTFin1}) in Example \ref{ex1} for  $h=0.02$ and $k=1,2,\ldots,5$}\label{fig1}	
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[width=10.cm]{fig2.eps}
%\includegraphics[width=10.cm]{figure2.eps}
\medskip
\caption{\rm\em Relative errors for methods (\ref{RecJTFin1}) in Example \ref{ex1} for $h=0.01$ and $k=1,2,\ldots,5$}\label{fig1a}	
\end{figure}

We consider now the actual errors, $|y_{n+k}-y(x_{n+k})|$ at a fixed point $x=x_{n+k}$  obtained by using the $k$-step method  (\ref{RecJTFin1}) for different stepsize $h$ and different $k$. We take $x=-0.5$ and $h=0.05$, $0.02$ and $0.01$. The corresponding errors are presented in Table~\ref{table1}. Numbers in parentheses indicate  decimal exponents.
\begin{table}[!htb]
\centering
\begin{tabular}{|c||c|c|c|c|c|}
\hline
$h$ & $k=1$ & $k=2$ & $k=3$ & $k=4$ & $k=5$ \\
\hline
$0.05$ & $2.27(-1)$ & $1.14(-1)$ & $4.76(-3)$  & $6.68(-5)$  & $1.45(-5)$ \\
$0.02$ & $4.10(-1)$ & $7.87(-2)$ & $1.20(-3)$ & $1.18(-5)$ & $6.89(-7)$   \\
$0.01$ & $5.46(-1)$ & $5.77(-2)$ & $4.21(-4)$ & $2.57(-6)$ & $6.34(-8)$ \\
\hline
\end{tabular}
\medskip
\caption{\rm\em Absolute  errors in the obtained sequences
$\{y_{n+k}\}_n$ at the point $x=x_{n+k}=-0.5$, using the $k$-step method (\ref{RecJTFin1}) for $k=1,2,\ldots,5$ and $h=0.05$, $0.02$ and $0.01$}\label{table1}\end{table}

Note that the effect of reducing the step-size $h$ to the accuracy of the solution is greater if $k$ is higher. Assuming an asymptotic relation in the form
\begin{equation}\label{AcEr}
e(h,k,x)=|y_{n+k}^h -y(x_{n+k})| \approx C_k h^{\alpha_k},	
\end{equation}
where $x_{n+k}=-1+(n+k)h=x=\mbox{const}$, and $C_k$ and $\alpha_k$ are some constants, we can calculate
the following quotient for two different steps $h_1$ and $h_2$,
\[\frac{e(h_1,k,x)}{e(h_2,k,x)}\approx  \left(\frac{h_1}{h_2}\right)^{\alpha_k}.\]
Therefore,
\begin{equation}\label{parALF}
\alpha_k=\frac{\log(e(h_1,k,x)/e(h_2,k,x)}{\log(h_1/h_2)}.	
\end{equation}

These values are  presented in Table~\ref{table1a} for $h_1/h_2=2$ and $h_1/h_2=5$. As we can see, the obtained values for the exponents $\alpha_k$ are very similar in these two cases.
\begin{table}[!htb]
\centering
\begin{tabular}{|c||c|c|c|c|c|}
\hline
$h_1/h_2$ & $\alpha_1$ & $\alpha_2$ & $\alpha_3$ & $\alpha_4$ & $\alpha_5$\\
\hline
$2$ & $-0.41$ & $0.45$ & $1.51$ & $2.20$ & $3.44$ \\
$5$ & $-0.54$ & $0.42$ & $1.51$ & $2.02$ & $3.38$\\
\hline
\end{tabular}
\medskip
\caption{\rm\em The parameters $\alpha_k$ obtained from (\ref{parALF}) for $k=1,2,\ldots,5$ and different stepsizes }\label{table1a}\end{table}


As we can see, the actual order of the method is reduced approximately for one order of magnitude. This effect is mentioned in Section~\ref{sec3} after Proposition~\ref{Prop1}.}
\end{example}


\begin{example} {\rm Here we consider again the Cauchy problem of Jacobi type
\[(1 - x^2)y' -2 xy =1 - x - 4 x^2 - 5 x^3+x y,\quad y(- 1) = 1,\]
whose exact solution is $y(x)=x^2 + x + 1$. According to Proposition~\ref{Prop1}, in this simple  case, we have
\[g(x)=G(x,y(x))=-4 x^3-3 x^2+1,\]
as well as
\[g'(x)=-12 x^2-6 x,\quad g''(x)=-24 x-6,\quad g'''(x)=-24,\quad g^{(iv)}(x)=0. \]

Now, we apply the $k$-step method \eqref{RecJTFin1} for $k=1(1)4$, $h=0.05$ and $h=0.01$ (see Tables \ref{ex1a} and \ref{ex1b}), and for starting values (when $k>1$) we use the exact values of the solution. In these tables, we only give the relative errors
of the obtained values for $x=-0.9,-0.8,-0.7,-0.6,-0.5$ (m.p. is machine precision).
\begin{table}[!htb]
\centering
\begin{tabular}{|c|c|c|c|c|}
\hline
$x$ & $k=1$ & $k=2$ & $k=3$ & $k=4$ \\
\hline
$-0.9$ & $5.62(-2)$ & $5.09(-3)$ &  &   \\
$-0.8$ & $4.88(-2)$ & $6.37(-3)$ & $3.17(-4)$ & m.p.  \\
$-0.7$ & $4.37(-2)$ & $7.03(-3)$ & $4.46(-4)$ & m.p. \\
$-0.6$ & $3.80(-2)$ & $7.29(-3)$ & $5.62(-4)$ & m.p. \\
$-0.5$ & $3.18(-2)$ & $7.23(-3)$ & $6.42(-4)$ & m.p.  \\
%$-0.4$ & $2.55(-2)$ & $6.90(-3)$ & $7.06(-4)$ & m.p. \\
%$-0.3$ & $1.96(-2)$ & $6.33(-3)$ & $7.56(-4)$ & m.p. \\
%$-0.2$ & $1.46(-2)$ & $5.58(-3)$ & $7.94(-4)$ & m.p. \\
%$-0.1$ & $1.11(-2)$ & $4.73(-3)$ & $8.23(-4)$ & m.p.  \\
\hline
\end{tabular}
\medskip
\caption{\rm\em Relative errors in the obtained approximate sequences $\{y_{n+k}\}_n$ using $k$-step methods \eqref{RecJTFin1} for $k=1(1)4$ and $h=0.05$}\label{ex1a}
\end{table}

\begin{table}[!htb]
\centering

\begin{tabular}{|c|c|c|c|c|}
\hline
$x$ & $k=1$ & $k=2$ & $k=3$ & $k=4$ \\
\hline
$-0.9$ & $9.81(-3)$ & $2.61(-4)$ & $3.19(-6)$ & m.p.  \\
$-0.8$ & $9.19(-3)$ & $2.82(-4)$ & $3.92(-6)$ & m.p. \\
$-0.7$ & $8.37(-3)$ & $2.94(-4)$ & $4.53(-6)$ & m.p. \\
$-0.6$ & $7.33(-3)$ & $2.96(-4)$ & $5.11(-6)$ & m.p. \\
$-0.5$ & $6.15(-3)$ & $2.89(-4)$ & $5.62(-6)$ & m.p.  \\
%$-0.4$ & $4.92(-3)$ & $2.72(-4)$ & $6.05(-6)$ & m.p. \\
%$-0.3$ & $3.77(-3)$ & $2.46(-4)$ & $6.39(-6)$ & m.p. \\
%$-0.2$ & $2.81(-3)$ & $2.14(-4)$ & $6.65(-6)$ & m.p. \\
%$-0.1$ & $2.14(-3)$ & $1.78(-4)$ & $6.85(-6)$ & m.p.  \\
\hline
\end{tabular}
\medskip
\caption{\rm\em Relative errors in the obtained approximate sequences $\{y_{n+k}\}_n$ using $k$-step methods \eqref{RecJTFin1} for $k=1(1)4$ and $h=0.01$}\label{ex1b}
\end{table}

As in Example \ref{ex1}, we consider  asymptotic relation (\ref{AcEr}) at $x=x_{n+k}=-0.5$, when $k=1,2,3$ and $h_1=0.05$ and $h_2=0.01$. The values of the corresponding exponent (\ref{parALF}) are presented in Table~\ref{tabb}.
\begin{table}[!htb]
\centering
\begin{tabular}{|c||c|c|c|c|c|}
\hline
$k$ & $k=1$ & $k=2$ & $k=3$ \\
\hline
$h_1/h_2=5$ & $\alpha_1=1.00$ & $\alpha_2=2.00$ & $\alpha_3=2.94$ \\
\hline
\end{tabular}
\medskip
\caption{\rm\em The parameters $\alpha_k$ obtained from (\ref{parALF}) for $k=1,2,3$ and $h_1/h_2=5$ }\label{tabb}\end{table}

As we can see, in this polynomial case, there is not   previously mentioned defect in the order. Note that  the local truncation error (\ref{divdf}) is equal to zero for each $k\ge 4$, because of $g^{(k)}(x)=0$. Also, we see that the actual errors in Tables \ref{ex1a} and \ref{ex1b} for $k=4$ are on the level of machine precision.
}
\end{example}



\begin{example} {\rm Now we consider the Cauchy problem of Laguerre type
\[xy'+(1-x)y=\frac{3x^2+1}{(x^2+1)^2}{\E}^{-x}y^2,\quad y(0)=1,\]
whose exact solution is $y=(x^2+1){\E}^x$. We apply the $k$-step method \eqref{RecJTFin} for $k=1(1)6$.   The corresponding relative errors
for  $h=0.05$ are given in Table \ref{ex31}, and for $h=0.01$ in Table \ref{ex32}. In these tables, we only give   relative errors of the obtained values for $x=0(0.1)1$. As in Example~\ref{ex1}, for starting values we use the exact values of   solution.
\begin{table}[!htb]
\centering
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
$x$ & $k=1$ & $k=2$ & $k=3$ & $k=4$ & $k=5$ & $k=6$\\
\hline
$0.1$ & $5.66(-2)$ & $3.74(-3)$ &  &  &  & \\
$0.2$ & $1.18(-1)$ & $1.43(-2)$ & $7.17(-4)$ & $1.89(-5)$ &   & \\
$0.3$ & $1.81(-1)$ & $2.62(-2)$ & $1.62(-3)$ & $7.55(-5)$ & $2.89(-6)$ & $6.33(-8)$ \\
$0.4$ & $2.45(-1)$ & $3.94(-2)$ & $2.62(-3)$ & $1.32(-4)$ & $5.83(-6)$ & $2.39(-7)$ \\
$0.5$ & $3.08(-1)$ & $5.42(-2)$ & $3.76(-3)$ & $1.98(-4)$ & $9.39(-6)$ & $4.16(-7)$ \\
$0.6$ & $3.70(-1)$ & $7.10(-2)$ & $5.10(-3)$ & $2.74(-4)$ & $1.33(-5)$ & $6.11(-7)$\\
$0.7$ & $4.30(-1)$ & $9.01(-2)$ & $6.66(-3)$ & $3.63(-4)$ & $1.79(-5)$ & $8.48(-7)$\\
$0.8$ & $4.88(-1)$ & $1.11(-1)$ & $8.49(-3)$ & $4.67(-4)$ & $2.33(-5)$ & $1.11(-6)$ \\
$0.9$ & $5.42(-1)$ & $1.35(-1)$ & $1.06(-2)$ & $5.88(-4)$ & $2.95(-5)$ & $1.42(-6)$\\
$1.0$ & $5.93(-1)$ & $1.62(-1)$ & $1.31(-2)$ & $7.29(-4)$ & $3.67(-5)$ & $1.77(-6)$\\
\hline
\end{tabular}
\medskip
\caption{\rm\em Relative errors in the obtained  sequences $\{y_{n+k}\}_n$ using $k$-step methods \eqref{RecJTFin} for $k=1,2,\ldots,6$ and $h=0.05$}\label{ex31}
\end{table}

\begin{table}[!htb]
\centering
\begin{tabular}{|c|c|c|c|c|c|l|}
\hline
$x$ & $k=1$ & $k=2$ & $k=3$ & $k=4$ & $k=5$ & \ \ $k=6$\\
\hline
$0.1$ & $5.26(-2)$ & $1.86(-3)$ & $2.63(-5)$ & $2.92(-7)$ & $2.97(-9)$ & $2.85(-11)$ \\
$0.2$ & $1.05(-1)$ & $4.17(-3)$ & $6.20(-5)$ & $7.37(-7)$ & $8.17(-9)$ & $8.70(-11)$ \\
$0.3$ & $1.58(-1)$ & $6.73(-3)$ & $1.02(-4)$ & $1.23(-6)$ & $1.39(-8)$ & $1.51(-10)$  \\
$0.4$ & $2.12(-1)$ & $9.66(-3)$ & $1.47(-4)$ & $1.80(-6)$ & $2.04(-8)$ & $2.24(-10)$ \\
$0.5$ & $2.67(-1)$ & $1.31(-2)$ & $2.00(-4)$ & $2.45(-6)$ & $2.80(-8)$ & $3.09(-10)$  \\
$0.6$ & $3.23(-1)$ & $1.71(-2)$ & $2.63(-4)$ & $3.23(-6)$ & $3.70(-8)$ & $4.10(-10)$ \\
$0.7$ & $3.80(-1)$ & $2.18(-2)$ & $3.38(-4)$ & $4.16(-6)$ & $4.76(-8)$ & $5.28(-10)$\\
$0.8$ & $4.35(-1)$ & $2.73(-2)$ & $4.26(-4)$ & $5.25(-6)$ & $6.02(-8)$  & $6.67(-10)$ \\
$0.9$ & $4.89(-1)$ & $3.38(-2)$ & $5.30(-4)$ & $6.53(-6)$ & $7.49(-8)$ & $8.32(-10)$\\
$1.0$ & $5.41(-1)$ & $4.12(-2)$ & $6.51(-4)$ & $8.03(-6)$ & $9.22(-8)$ & $1.02(-9)$ \\
\hline
\end{tabular}
\medskip
\caption{\rm\em Relative errors in the obtained  sequences $\{y_{n+k}\}_n$ using $k$-step methods \eqref{RecJTFin} for $k=1,2,\ldots,6$ and $h=0.01$}\label{ex32}
\end{table}

Using Proposition~\ref{Prop1} we determine the
 main term of the truncation error at the point
$x_{n+k}=x$, for example, when $h=0.01$ and $x=0.5$.

Since
\[\int_0^1\binom{\lambda+k-1}{k}{\E}^{-(x-h+\lambda h)}{\D}\lambda
={\E}^{-x}Q_k,\quad Q_k=\frac{1}{k!}\int_0^1 (\lambda)_k {\E}^{-(\lambda-1)h}{\D}\lambda,\]
we first calculate  the   values:
$Q_1=0.501671$, $Q_2=0.41792$, $Q_3=0.376058$, $Q_4=0.34955$, $Q_5=0.330639$, and  $Q_6=0.316389$.
The corresponding derivatives %$g^{(k)}(x)$ are
\[g^{(k)}(x)=x y^{(k+1)}+(k+1-x)y^{(k)}-ky^{(k-1)}\]
are
\begin{align*}
 g'(x)&={\E}^x \left(3 x^2+6 x+1\right),&	\qquad \ \ \, g''(x)&={\E}^x \left(3 x^2+12 x+7\right),\\
 g'''(x)&={\E}^x\left(3 x^2+18x+19\right),&\quad \
g^{(iv)}(x)&={\E}^x \left(3 x^2+24 x+37\right),\\
 g^{(v)}(x)&={\E}^x \left(3 x^2+30 x+61\right),&\quad g^{(vi)}(x)&={\E}^x\left(3x^2+36 x+91\right),
\end{align*}
where $\xi_k\in \bigl(x-(k+1)h,x-h\bigr)$, $k=1,2,\ldots,6$.

Now, taking $x=x_{n+k}=(n+k)h=0.5$, $h=0.01$, and $\xi_k=x-h=0.49$ in (\ref{divdf}), we obtain an approximation  of the main term $(\widehat T_h^w)_{n+k}$ in the form
\[(\widehat T_h^w)_{n+k} \approx \frac{h^k g^{(k)}(x-h)Q_k}{x},\qquad k=1,2,\ldots,6,\]
whose numerical values for $k=1,2,\ldots,6$, after dividing by $y(0.5)=2.0609$, are: $3.70(-2)$, $9.00(-4)$, $1.70(-5)$, $2.74(-7)$, $4.00(-9)$, $5.48(-11)$, respectively. As   expected,
the actual global errors from Table \ref{ex32} (the row referring to $x=0.5$)   are larger compared to  the corresponding local truncation errors.


\begin{table}[!htb]
\centering
\begin{tabular}{|c||c|c|c|c|c|c|}
\hline
$x$ & $\alpha_1$ & $\alpha_2$ & $\alpha_3$ & $\alpha_4$ & $\alpha_5$ & $\alpha_6$\\
\hline
$0.5$ & $0.08$ & $0.95$ & $1.92$ & $2.87$ & $3.81$ & $4.75$ \\
$1.0$ & $0.05$ & $0.93$ & $1.94$ & $2.90$ & $3.86$ & $4.80$\\
\hline
\end{tabular}
\medskip
\caption{\rm\em The parameters $\alpha_k$ obtained from (\ref{parALF}) for $k=1,2,\ldots,6$ and $h1/h2=2$ at two points $x=0.5$ and $x=1.0$.}\label{tabela3}
\end{table}

Finally, as in Example~\ref{ex1}, we assume the behavior of the  actual errors in the form (\ref{AcEr}), where  $x_{n+k}=(n+k)h=x=\mbox{const}$. We compare actual errors obtained for $h_1=0.02$ and $h_2=0.01$ at two points $x=0.5$ and $x=1.0$. The results for $\alpha_k$, $k=0,1,\ldots,6$ are presented in Table~\ref{tabela3}.  As we can see, the obtained values of the exponents $\alpha_k$ at these points are very close, but again with a defect of one order in its magnitude!
}
\end{example}


%%%% Acknowledgment %%%%%%%%
\section*{Acknowledgements}
The authors are deeply grateful to the anonymous referee for his/her
 comments and constructive suggestions for improvements of
 this paper.

 The work of M. Masjed-Jamei has been supported by the Alexander von Humboldt Foundation under the grant number:
 Ref 3.4--IRN--1128637--GF--E.

 The work of G.V. Milovanovi\'c  was supported in part by
  the Serbian Academy of Sciences and Arts (No.\ $\Phi$-96) and by the
  Serbian Ministry of Education, Science and Technological Development
  (No.\ \#OI 174015).


%%%% Bibliography  %%%%%%%%%%
\begin{thebibliography}{99}

\bibitem{3}
 {\sc K.\,E.\,Atkinson, W.\,Ham, D.\,Stewart}, {\em Numerical Solution of Ordinary Differential Equations},  John Wiley \& Sons Inc., Hoboken, NJ, 2009.

\bibitem{Bash}
{\sc F.\,Bashforth, J.\,C.\,Adams},  {\em An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid, with an explanation of the method of integration employed in constructing the tables which give the theoretical forms of such drops}, Cambridge University Press, Cambridge, 1883.


\bibitem{Butcher2000}
 {\sc J.\,C.\,Butcher}, {\em Numerical methods for ordinary differential equations in the 20th century}, J. Comput. Appl. Math. {\bf125}(2000), 1--29.

\bibitem{6}
 {\sc J.\,C.\,Butcher}, {\em  Numerical Methods for Ordinary Differential Equations}, John Wiley \& Sons, Ltd, Chichester, 2016.

 \bibitem{7}
 {\sc P.\,L.\,Butzer, M.\,Hauss}, {\em On Stirling functions of the second kind}, Stud. Appl. Math.  {\bf 84}(1991),   71--91.

 \bibitem{8}
 {\sc G.\,Dahlquist}, {\em  Convergence and stability in the numerical integration of ordinary differential equations},  Math. Scand.
 {\bf4}(1956), 33--53.

 \bibitem{WG}
 {\sc W.\,Gautschi}, {\em Numerical Analysis}, 2nd ed., Birkh\"auser,
 New York, 2012.

\bibitem{20}
 {\sc M.\,Masjed-Jamei},  {\em On constructing  new interpolation formulas using linear operators and an operator type of quadrature rules},  J. Comput.  Appl. Math. {\bf 216}(2008),  307--318.

 \bibitem{22}
 {\sc  M.\,Masjed-Jamei, G.V. Milovanovi\'c, M.A. Jafari}, {\em Explicit forms of weighted quadrature rules with geometric nodes},   Math. Comput. Modelling {\bf 53}(2011), 1133--1139.

  \bibitem{Mas_Mil2008}
{\sc G.\,Mastroianni, G.\,V.\,Milovanovi\'c}, {\em Interpolation Processes -- Basic Theory and Applications}, Springer-Verlag, Berlin-Heidelberg, 2008.

\bibitem{gvm88}
{\sc G.\,V.\,Milovanovi\'c}, {\em Numerical Analysis III}, Nau\v cna Knjiga, Belgrade, 1988. (in Serbian)

  \bibitem{Moulton}
{\sc F.\,R.\,Moulton}, {\em  New Methods in Exterior Ballistics}, University of Chicago, Chicago, 1926.

  \bibitem{18}
{\sc J.\,Sandor, B.\,Crstici}, {\em Handbook of Number Theory}, Part II, Springer/Kluwer Academic Publishers, Dordrecht, 2004.



 \end{thebibliography}
%\end{linenumbers}
\end{document}