\documentclass[12pt,twoside]{article}
\usepackage{latexsym,amsmath,amssymb}
\renewcommand{\baselinestretch}{1}
\usepackage{graphics}
\setcounter{page}{1} \setcounter{page}{1}
\setlength{\textheight}{21.6cm}

\setlength{\textwidth}{14cm}

\setlength{\oddsidemargin}{1cm}

\setlength{\evensidemargin}{1cm}

\pagestyle{myheadings}

\thispagestyle{empty}
  \markboth{\small{  Parviz Darania  and Fatemeh Sotoudehmaram }}
  {\small{On the linear stability and super convergence
of  multi step...
 }}

\date{}

\begin{document}

%\centerline{\bf Applied Mathematical Sciences, Vol. x, 200x,
%no.xx, xxx - xxx}

%\centerline{}

%\centerline{}



\centerline {\Large{\bf   On the linear stability and super convergence
of
  }}
\centerline{}

\centerline{\Large{\bf  multi step collocation solutions  for    functional
 }}
 \centerline{}

\centerline{\Large{\bf Volterra integral equations  }}


\centerline{} \centerline{} \centerline{\bf {Parviz  Darania$^1$ and Fatemeh
Sotoudehmaram$^2$}}

 \centerline{}



 \centerline{}

 \centerline{  $^{1,2}$Department of Mathematics,
Faculty of Science, } \centerline{Urmia University, P.O.Box 165,
Urmia-Iran}

 \centerline{}
\centerline{ $^1$ p.darania@urmia.ac.ir}
\centerline{$^2$fatemeh\_sotoudehmaram@yahoo.com,}
  \centerline{}
 \centerline{\bf {}}


%################################################
\newtheorem{defi}{Definition}[section]
\newtheorem{exam}[defi]{Example}
%\theoremstyle{plain}
\newcommand{\dst}{\displaystyle}
\newtheorem{lemma}[defi]{Lemma}
\newtheorem{pro}[defi]{Proposition}
\newtheorem{theorem}[defi]{Theorem}
\newtheorem{cor}[defi]{Corollary}
%\theoremstyle{remark}
\newtheorem{rem}[defi]{Remark}
%\newtheorem*{remn}{Remark}
\newtheorem{case}{Case:}
\newcommand{\proof}{\textbf{Proof : }}
\newcommand{\pa}{$\square$}
\newcommand{\pp}{\square}
\newcommand{\ka}{\kappa_a}
\newcommand{\bdm}{\begin{displaymath}}
\newcommand{\edm}{\end{displaymath}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\f}[2]{\frac{#1}{#2}}
\newcommand{\si}{\sigma}
\newcommand{\co}{\nabla^\xi}
\newcommand{\Co}{\nabla}
\newcommand{\rd}{$\Bbb{R}^2$}
\newcommand{\rc}{$\Bbb{R}^3$}
\newcommand{\rn}{$\Bbb{R}^n$}
\newcommand{\rr}{\Bbb{R}}
\newcommand{\n}{\Bbb{N}}
%\newcommand{\c}{\Bbb{C}}
\newcommand{\p}{\Bbb{\Pi}}
\newcommand{\ba}{\bar{M}}
\newcommand{\all}{\forall}
\newcommand{\no}{\nonumber}
\newcommand{\ben}{\begin{enumerate}}
\newcommand{\een}{\end{enumerate}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract}
This paper expands on a topic that was explored in recent research by P. Darania et al. (2020), in which the multi-step collocation method (MSCM) has been studied for solving FIEs along with their convergence analysis. This study focuses on the analysis of the linear stability and super-convergence of MSCM for a general form of delay function  $\Theta(t)$, which includes discontinuity points where the solutions exhibit low regularity. In addition, we present the stability region and the super-convergence order with various examples, thereby demonstrating the accuracy and efficiency of this method. Moreover, the super-convergence and stability regions of the MSCM for FIEs with delay, $\Theta (t)=t-\tau$, that do not exhibit any points of discontinuity have been investigated in \cite{darania-pishbin}.

\end{abstract}

{\bf   2010 Mathematics Subject Classification:}  45G05; 45D05; 65L05; 65L20.\\

\vspace{0.5cm} {\it \bf  Keywords:}  Functional Volterra
 integral equations, Multi-step collocation method,
Super-convergence, Linear stability analysis.
% ----------------------------------------------------------------------
\goodbreak
\section{Introduction}

In the description of certain phenomena in physics, it is required to take into account that the rate of propagation of interaction is finite, such as relativistic dynamics and nuclear reactors. Furthermore, it should be noted that these phenomena can usually be modeled with functional differential equations, which can be easily converted to FIEs by integrating them. This conversion can sometimes be useful from a numerical point of view due to the advantages of the numerical stability properties of FIEs when compared to functional differential equations \cite{lin}.


 Let us consider the   FIEs of the form

\begin{equation}
  \label{2*}
  y(t)=\left \{
  \begin{array}{lll}
   \phi(z)+\dst\int_{z_{0}}^{t}k_{1}(z,\tau,y(\tau))d\tau+\dst \int_{z_{0}}^{ \Theta(z)}
k_{2}(z,\tau,y(\tau))d\tau,& z\in (z_{0},Z],\\
    &\\
   g(z),&t\in[\theta(z_0),z_0],
    \end{array}\right.
   \end{equation}

Here, $k_{1} \in C(\Omega)$,  $\Omega=\lbrace (z,\tau) : z_{0}\leq \tau \leq z \leq
Z,  z \in J :=[z_{0},Z] \rbrace$ and $k_{2}$ is supposed to be
continuous in  $\Omega_{\Theta} =\lbrace(z,\tau) : \theta(z_{0}) \leq \tau \leq
\Theta(z), z \in J \rbrace$ and $g,\phi, $ are at least continuous on
their domains.

The following conditions will apply to the delay function $\Theta$ within this paper:
\begin{itemize}
\item[(A1)] $\Theta(z)=z-\varphi(z)$, with $\varphi(z) \geq \varphi_{0}>0$ for $z \in
J$,
\item[(A2)] $\Theta$ is strictly increasing on $J$,
\item[(A3)] $\varphi\in C^{d}(J)$ for some $d \geq 0$.
\end{itemize}

The non-vanishing delay $\Theta(z)$  induces the primary discontinuity points $\{\xi_{\mu}\} $ for the solution of
\eqref{2*}: they are ascertained by the recursion
$$\Theta(\xi_{\mu})=\xi_{\mu-1}, \qquad \mu \geq 0, \quad \xi_{-1}=\Theta(z_{0}),~\xi_{0}=z_{0}.$$

These conditions ensure that the given discontinuity points possess uniform separation property
$$\xi_{\mu}-\xi_{\mu-1}=\varphi(\xi_{\mu}) \geq \varphi_{0} > 0, \qquad \mbox{for all} \quad \mu \geq 0.$$

\begin{theorem}\cite{bookbrunner}
\label{th1}   Assume that the provided functions in the equation
\begin{equation}
  \label{2-1*}
  y(z)=\left \{
  \begin{array}{lll}
    \phi(z)+\dst\int_{z_{0}}^{z}k_{1}(z,\tau)y(\tau)d\tau+\dst \int_{z_{0}}^{ \Theta(z)}
k_{2}(z,\tau)y(\tau)d\tau,& z\in (z_{0},Z],\\
    &\\
    g (z),&z\in[\Theta(z_0),z_0],
    \end{array}\right.
   \end{equation}
exhibit continuity within their respective domains and the delay function $\Theta$ holds the above conditions  $(A1)-(A3)$.
 Therefore, there exists a unique and bounded $y \in C(z_0,\Omega ]$ for any initial function  $g \in C[\Theta (z_0),z_0]$ which solves the delay integral equation \eqref{2-1*} on the interval $(z_0,Z]$  and corresponds to $ g$ on the interval $[\Theta (z_0),z_0]$.
 In general, this solution possesses a finite (jump) discontinuity at $z = z_0$:

$$
 \lim _{z\rightarrow z^+_0}y(z)\neq \lim_{z\rightarrow z^-_0} y(z) =g(t_0).
$$
The solution will be continuous at $z = z_0$ if, and only if, the initial function satisfies the following equation:
$$
\phi(z_0) - \int_{\Theta(z_0)}^{ z_0} k_2(z_0, \tau)g(\tau)d\tau = g(t_0).
$$
\end{theorem}

Due to the fact that solutions of delay problems with non-vanishing delays generally experience loss of regularity at the primary discontinuity points $\{ \xi_{\mu} \}$,
the mesh $J_{h}$ underlying the collocation space must incorporate these points if the collocation solution is to achieve its optimal global (or local) order (of super-convergence). Thus, for
$J:=(z_{0},Z]$ and for $\mu=0, \ldots ,M$, we will utilize meshes of the form

\begin{equation}\label{5*}
J_{h}: = \bigcup _{\mu=0}^M J_{h}^{(\mu)}, \qquad
J_{h}^{(\mu)}:=\{z_{n}^{(\mu)}~:~
\xi_{\mu}=z_{0}^{(\mu)}<z_{1}^{(\mu)}< \cdots
<z_{N}^{(\mu)}=\xi_{\mu+1}\},
\end{equation}
and for $\mu=-1$,
$$J_{h}^{(-1)}:=\{z_{n}^{(-1)}~ : ~ \Theta(z_{0})=z_{0}^{(-1)} <z_{1}^{(-1)} < \cdots < z_{N}^{(-1)}=z_{0} \}.$$

\begin{defi}
A mesh $J_{h}$ for $J:=[z_{0},Z]$ is considered to be $\Theta$-invariant
if it is constrained, and if
\begin{equation}\label{6*}
\Theta(J_{h}^{(\mu)})=J_{h}^{(\mu -1)}, \qquad \qquad \mu =1, \ldots
, M,
\end{equation}
holds.
\end{defi}
 So, we will have the following rigorous notations:
 \begin{equation*}
  \sigma_{n}^{(\mu)}:=(z_{n}^{(\mu)},z_{n+1}^{(\mu)}],\quad \bar{\sigma}_{n}:=[z_{n}^{(\mu)},z_{n+1}^{(\mu)}],\quad h_{n}^{(\mu)}:=z_{n+1}^{(\mu)}-z_{n}^{(\mu)},\\
 \end{equation*}
 \begin{equation*}
 h^{(\mu)}:=\max\limits_{n}h_{n}^{(\mu)}, \qquad h:=\max\limits_{(\mu)}h^{(\mu)},\qquad 0\leq n \leq N-1.
 \end{equation*}
Considering a $\Theta$-variant mesh $J_{h}$, the collocation solution $v_{h}$ will be a component of a piecewise polynomial space

\begin{equation}\label{7*}
S_{m+d}^{(d)}(J_{h}):=\{u \in C^{d}(J_{h}) : u|_{\sigma_{n}^{(\mu)}}
\in \Pi_{m+d}, 0 \leq n \leq N-1,~0 \leq \mu \leq M\}.
\end{equation}
Let's choose the collocation points as
\begin{equation}\label{8*}
X_{h}:= \bigcup_{\mu=0}^{M} X_{h}^{(\mu)},
\end{equation}
These points are determined based on the $M+1$ local sets

$$X_{h}^{(\mu)}=\{z_{n,i}^{(\mu)}=z_{n}^{(\mu)}+c_{i}h_{n}^{(\mu)} : 0<c_{1}< \cdots <c_{m} \leq 1, ~n=0,..., N-1\}.$$
It’s obvious that the cardinality of $X_{h}$ will be equivalent to $Nm$. For more detail, see \cite{bookbrunner} and references therein.

The representation of the collocation solution $v_{h}\in S_{m-1}^{(-1)}(J_{h})$ for the delay  IEs \eqref{2*} is given by the collocation equation:

\begin{equation}\label{16*}
 v_{h}(t)=\left \{
\begin{array}{lr}
\phi(z)+\dst\int_{z_{0}}^{z}k_{1}(z,\tau,v_{h}(\tau))d\tau+\dst \int_{z_{0}}^{
\Theta(z)}
k_{2}(z,\tau,v_{h}(\tau))d\tau, & z \in X_{h},\\
g(z), &z \in [\Theta(z_{0}),z_{0}],
\end{array}\right.
\end{equation}
on $\sigma_{n}^{(\mu)}$. The lag function $\Theta=\Theta(z)=z-\varphi(t)$ will be assumed to satisfy the conditions $ (A1)-(A3)$. Additionally, we will consider the mesh  $J_{h}$ on the interval $J:=[z_{0},Z]$. The MSCM is derived from incorporating $r$ previous time steps into the collocation polynomials. Particularly, we will seek a collocation polynomial with a restriction to the interval $\sigma_{n}^{(\mu)}$. Proofs of most of the results can be found in \cite{darania-sutudeh, darania-mara,darania-tabriz, bookbrunner}. In \cite{bookbrunner}, Brunner employed collocation-type methods to compute numerical solutions for nonlinear delay IEs and investigated their association with iterated collocation methods.
In \cite{darania-mara} and \cite{darania-tabriz}, in order to achieve good convergence, local super-convergence, and stability properties, a class of MSCM was constructed by loosening some collocation conditions. The analysis of the MSCM for solving delay IEs with a delay function $t-\tau$ was first studied by P. Darania \cite{darania-sutudeh,darania-tabriz} and \cite{baker, brunner1}. Moreover, a general class of MSCM that relies on $r$ fixed number of previous time steps and $m$ collocation points was defined by him in \cite{darania-mara}.
These methods exhibit a uniform order of $m+r$ for all choices of collocation parameters, with a local super-convergence order of $2m+2r-1$ in mesh points for a specific selection of collocation parameters. In contrast, classical collocation methods possess a uniform order of $m$ for all choices of collocation parameters and local super-convergence properties at the mesh points, with an order of $2m-2$ for Gauss and Lobatto points and an order of $2m-1$ for Radau II points (see \cite{bookbrunner}). The construction of a MSCM to solve a general class of nonlinear delay IEs, including two types of linear and nonlinear lag functions $
\Theta(t)$, and the investigation of the convergence analysis of this method were carried out by P. Darania and F. Shotoudehmaram in \cite{darania-sutudeh}. In this paper, our focus is on investigating the super-convergence with an extensive stability region of the MSCM for solving delay integral equations. The next sections of this paper are arranged as follows: In Section 2, the orders of super-convergence and linear stability of this technique are ascertained, and the manuscript is concluded in Section 3, by demonstrating the effectiveness of this technique through some numerical instances.

%\section{ Main Results}
 \section{ Super-convergence}
Here, we present the formula for local super-convergence of equation  \eqref{16*} as well as the implementation method in the MSCM.
In order to prove the super-convergence of the solution to the problem \eqref{16*}, the super-convergence analysis will be discussed using the method proposed in \cite{darania-sutudeh}. We will also refer to the theorem that is proved in \cite{darania-sutudeh} and represents the convergence result of the problem \eqref{16*}.


\begin{theorem}\cite{darania-sutudeh}
\label{th2}
  Let us consider the function in \eqref{2-1*} satisfying $\phi \in C^{m+r}(I)$, $ k_{1} \in C^{m+r}(\Omega)$, $ k_{2} \in C^{m+r}(\Omega_{\Theta})$, $g \in C^{m+r}([z_{0},Z]$, and for $z \in [\Theta(z_{0}),z_{0}]$ the integral
\begin{equation}\label{57}
G(t):= \dst\int_{z_{0}}^{\Theta(z)}k_{2}(z,\tau)g(\tau)d\tau,
\end{equation}
is determined exactly. Also, suppose that the initial error for any $z\in (\xi_{\mu} ,\xi_{\mu+1} ],~\mu=0, \ldots, M$ is $\Vert y-v_{h}
\Vert_{\infty}=\mathcal{O}(h^{m+r}),$ and $ \rho(\mathbf{A})<1, $, where $\rho$ indicates the spectral radius and
\begin{equation}\label{59}
\pmb{A}=\begin{bmatrix}
\begin{array}{c|c}
\pmb{0}_{r-1,1} & \mathbf{I}_{r-1} \\
\hline \pmb{L_{r-1}(1)} & \pmb{L}_{r-2}(1), \ldots ,
\pmb{L}_{0}(1)
\end{array}
\end{bmatrix}.
\end{equation}
Then, for any sufficiently small $h$, the restricted mesh collocation solution \linebreak$ v_{h} \in S_{m-1}^{(-1)}(J_h) $
for \eqref{2-1*} satisfies $
\Vert \varepsilon \Vert_{\infty} \leq Ch^{m+r},$ where $
\varepsilon(t)=y(z) - v_h(z) $ represents the error of the exact collocation method and $C$ is a positive constant independent of $h$. This estimation holds for all collocation parameters  $ \{c_{i}\} $  with $ 0 \leq c_{1} \leq \cdots  \leq c_{m} \leq 1.$\\
\end{theorem}


\begin{theorem} \cite{darania-sutudeh}
\label{th3} Assuming that theorem \ref{th2} is satisfied, with the exception that the integrals

\begin{align*}
G(z)&= \dst\int_{z_{0}}^{\Theta(z)} k_{2}(z,\tau)g(\tau)d\tau,\qquad
z=z_{n,i}^{(\mu)},\quad  n=0,1, \ldots N-1,\quad \mu=0 ,
\end{align*}
are approximated by quadrature $\bar{G}(z)$, along with associated quadrature errors $E_{0}(t):=G(z) -\bar{G}(z)$, where $ \Vert E_{0}(z) \Vert \leq h^{q}, $ for some $q>0$.
In this case, the collocation solution  $v_{h} \in S_{m-1}^{(-1)}(J_h)$ satisfies for any sufficiently small $h>0$, where, $\Vert \varepsilon \Vert_{\infty}
\leq Ch^{p},$  with $ p=\min \{ m+r ,q\}$, and $C$ is a finite constant that is independent of $h$.

 \end{theorem}

Now, we will obtain the local super-convergence in the interior
points $J^{(\mu)}_{h}, ~\mu=0,1,\ldots, M$. For this purpose, in the
same manner in \cite{darania-sutudeh}, we define

\begin{equation}
\label{14*}
 v_{h}(z_{n}^{(\mu)}+\tau h_{n}^{(\mu)})= \dst \sum_{k=0}^{r-1}L_{k}(\tau)y_{n-k}^{(\mu )}+
 \sum_{j=1}^m \hat{L}_{j}(s)V_{n,j}^{(\mu)}, \qquad \tau \in (0,1] ,
\end{equation}
where $L_{k}(\tau)$ and $\hat{L}_{j}(\tau)$  are the following polynomials with a degree of $m + r - 1$ defined by

\begin{equation}\label{12*}
 \begin{array}{llll}
   L_k(\tau)=\displaystyle\prod_{i=1}^{m}\limits\frac{s-c_i}{-k-c_i}
    \cdot\displaystyle \prod _{\substack{i=0\\i\neq k}}^{r-1}\limits\frac{s+i}{-k+i},
    &
   \hat{L}(\tau)=\displaystyle\prod _{i=0}^{r-1}\limits\frac{s+i}{c_j+i} \cdot\displaystyle\prod _{\substack{i=1\\i\neq
    j}}^{m}\limits\frac{s-c_i}{c_j-c_i},
 \end{array}
\end{equation}
and the collocation parameters are supposed to satisfy  $c_{1}\neq 0$
and $c_{i} \neq c_{j}$ for $i\neq j$ and

\begin{equation}\label{15*}
V_{n,j}^{(\mu)}= v_{h}(z_{n,j}^{(\mu)}),\quad n=r, \ldots , N-1.
\end{equation}



  \begin{theorem}
\label{th4}
 Let the assumption of the theorem \ref{th2} holds with \linebreak $p=2m+r-1$ and let the collocation parameters $\{c_{i}\}_{i=1}^{m}, $ be the solution of the system
 \begin{equation}\label{107}\left\{
 \begin{array}{lr}
c_{m}=1,\\
\dfrac{1}{i+1}-
\dst\sum_{k=0}^{r-1}\beta_{k}(-k)^{i}-\dst\sum_{j=1}^{m}\gamma_{j}(c_{j})^{i}=0,
\quad i=m+r, \ldots , 2m+r-2,
 \end{array}\right.
 \end{equation}
 where
 \begin{align}\label{108}
 \beta_{k}=\dst\int_{0}^{1}L_{k}(\tau)d\tau, \qquad \gamma_{j}=\dst\int_{0}^{1}\hat{L}_{j}(\tau)d\tau.
 \end{align}
 Moreover, assume that the delay integral
 \begin{equation}\label{109}
G(z)=\dst\int_{z_{0}}^{\Theta(z)}k_{2}(z,\tau,g(\tau))d\tau,
 \end{equation}
 can be evaluated analytically. Then
 \begin{equation}\label{110}
 \max_{1\leq n \leq N}\vert \varepsilon(z_{n}) \vert \leq Ch^{p}.
 \end{equation}
 \end{theorem}
 \begin{proof}
 The collocation equation \eqref{2-1*} can be expressed in the form
 \begin{equation}\label{111}
 v_h(t)=-\delta_{h}(z)+\dst\int_{z_{0}}^{z}k_{1}(z,\tau)v_h(\tau)d\tau+\dst\int_{z_{0}}^{\Theta(z)}k_{2}(z,\tau)v_h(\tau)d\tau, \quad z \in J,
 \end{equation}
 where the defect function $\delta_h$ is zero at interior points $J^{(\mu)}_{h}$,
 \begin{equation}\label{112}
 \delta_h(t_{n,i}^{(\mu)})=0,
 \end{equation}
  Additionally, we get $\delta(z)=0$ for $z \leq z_{0}$. The collocation error $\varepsilon=y-v_h$ can be obtained by solving the integral equation
 \begin{align}\label{113}
  \varepsilon(z) &=\delta_{h}(z)+\dst\int_{z_{0}}^{z}k_{1}(z,\tau)\varepsilon_{h}(\tau)d\tau+F(z), \qquad z \in J,
 \end{align}
 where
 \begin{equation}\label{114}
 F(z)=\dst\int_{z_{0}}^{\Theta(z)}k_{2}(z,\tau)\varepsilon(\tau)d\tau, \qquad z \in [z_{0},Z].
 \end{equation}
For $z \in [\Theta(z_{0}),z_{0}]$, we obtain $F(z)=0$. Therefore, the error equation \eqref{113} simplifies to a classical Volterra equation, which has a unique solution given by
  \begin{equation}\label{115}
  \varepsilon(z)=\delta_h(z)+\dst\int_{z_{0}}^{z}R_{1}(z,\tau)\delta_h(\tau)d\tau,
  \end{equation}
where $R_{1}$ represents the resolvent kernel corresponding to the given kernel $k_{1}$. Since $Z=\xi_{M+1},$, where $M$ is a positive integer, and for $z \in [\xi_{\mu},\xi_{\mu+1}]$, $1 \leq \mu \leq M-1,$ the collocation error $\varepsilon(z)$ can be represented by equation \eqref{113} in the form
 \begin{equation}\label{116}
  \begin{array}{lr}
  \varepsilon(z) =\delta_{h}(z)+\dst\int_{\xi_{\mu}}^{z}R_{\mu}(z,\tau)\delta_{h}(\tau)d\tau
  +\dst\sum_{\kappa=0}^{\mu-1}\dst\int_{\xi_{\kappa}}^{\xi_{\kappa+1}} R_{\kappa}(z,\tau)\delta_{h}(\tau)d\tau\nonumber\\
\quad
+\displaystyle\int_{\xi_{\mu-1}}^{\Theta(z)}Q_{\mu-1}(z,\tau)\delta_{h}(\tau)d\tau
+\dst\sum_{\kappa=0}^{\mu-1}\dst\int_{\xi_{\kappa}}^{\xi_{\kappa+1}}Q_{v}(z,\tau)\delta_{h}(\tau)d\tau,\qquad
z \in J,
\end{array}
 \end{equation}
Here, $R_{\mu}(z,\tau)$ and $Q_{\kappa}(z,\tau)$ indicate functions that exhibit continuity on their respective domains and rely on $k_{1}$ and $k_{2}$.

 Now, for $z=z_{n}^{(\mu)}$, we obtain
\begin{equation}\label{117}
 \begin{array}{lr}
\varepsilon(z_{n}^{(\mu)})=\dst\sum_{l=0}^{n-1}h_{l}^{(\mu)}\dst\int_{0}^{1}R_{\mu}(z_{n}^{(\mu)}, z_{l}^{(\mu)}+\tau h_{l}^{(\mu)})\delta_{h}(z_{l}^{(\mu)}+\tau h_{l}^{(\mu)})d\tau\nonumber\\
 \quad +\dst\sum_{\kappa=0}^{\mu-1}\dst\sum_{l=0}^{N-1}h_{l}^{(\kappa)}\dst\int_{0}^{1}R_{\kappa}(z_{n}^{(\kappa)}, z_{l}^{(\kappa)}+\tau h_{l}^{(\kappa)})\delta_{h}(z_{l}^{(\kappa)}+\tau h_{l}^{(\kappa)})d\tau\nonumber\\
 \quad +\displaystyle{\sum_{l=0}^{n-1}h_{l}^{(\mu-1)}\int_{0}^{1}Q_{\mu-1}(z_{n}^{(\mu)},z_{l}^{(\mu-1)}+\tau h_{l}^{(\mu-1)})\delta_{h}(z_{l}^{(\mu-1)}+\tau h_{l}^{(\mu-1)})d\tau}\nonumber\\
\quad
+\dst\sum_{\kappa=0}^{\mu-1}\dst\sum_{l=0}^{N-1}h_{l}^{(\kappa)}\dst\int_{0}^{1}Q_{\kappa}(z_{n}^{(\kappa)},
z_{l}^{(\kappa)}+\tau h_{l}^{(\kappa)})\delta_{h}(z_{l}^{(\kappa)}+\tau h_{l}^{(\kappa)})d\tau.
\end{array}
\end{equation}

Consider the quadrature formula
\begin{equation}\label{118}
\dst\int_{0}^{1}f(\tau)d\tau \approx
\dst\sum_{k=0}^{r-1}\beta_{k}f(-k)+\dst\sum_{j=1}^{m}\gamma_{j}f(c_{j}),
\end{equation}
and assume that the collocation parameters  $\{c_{i}\}_{i=1}^{m}, $ are the solution of the system
 \begin{equation}\label{107}\left\{
 \begin{array}{lr}
c_{m}=1,\\
\dfrac{1}{i+1}-
\dst\sum_{k=0}^{r-1}\beta_{k}(-k)^{i}-\dst\sum_{j=1}^{m}\gamma_{j}(c_{j})^{i}=0,
\quad i=m+r, \ldots , p-1,
 \end{array}\right.
 \end{equation}
with
 \begin{align}\label{108}
 \beta_{k}=\dst\int_{0}^{1}L_{k}(\tau)d\tau, \qquad
 \gamma_{j}=\dst\int_{0}^{1}\hat{L}_{j}(\tau)d\tau.
 \end{align}
Let, $h^{(\mu)}:=\max_{n}\limits h_{n}^{(\mu)}$ and
$h:=\max_{\mu}\limits h^{(\mu)}$. Then, for the computation of the
integrals in \eqref{118}, we get
\begin{equation}\label{119}
 \begin{array}{rlr}
\varepsilon(z_{n}^{(\mu)})=&h\displaystyle\sum_{l=0}^{n-1}\Bigg{(}\displaystyle\sum_{k=0}^{r-1}\beta_{k}R_{\mu}(z_{n}^{(\mu)},z_{l}^{(\mu)}-kh_{l}^{(\mu)})\delta_{h}(zt_{l}^{(\mu)}-kh_{l}^{(\mu)})  \\

& +\displaystyle\sum_{j=1}^{m}\gamma_{j} R_{\mu}(z_{n}^{(\mu)},z_{l}^{(\mu)}+c_{j}h_{l}^{(\mu)})\delta_{h}(z_{l}^{(\mu)}+c_{j}h_{l}^{(\mu)})\Bigg{)} +h\displaystyle\sum_{l=0}^{n-1}E_{n,l,\mu}\\

&+h \displaystyle\sum_{\kappa=0}^{\mu-1} \displaystyle\sum_{l=0}^{N-1}\Bigg{(}\displaystyle\sum_{k=0}^{r-1}\beta_{k}R_{\kappa}(z_{n}^{(\mu)},z_{l}^{(\kappa)}-kh_{l}^{(\kappa)})\delta_{h}(z_{l}^{(\kappa)}-kh_{l}^{(\kappa)})\\
&+\displaystyle\sum_{j=1}^{m} \gamma_{j}R_{\kappa}(z_{n}^{(\mu)},z_{l}^{(\kappa)}+c_{j}h_{l}^{(\kappa)})\delta_{h}(z_{l}^{(\kappa)}+c_{j}h_{l}^{(\kappa)})\Bigg{)}+h\displaystyle\sum_{l=0}^{N-1}\displaystyle\sum_{\kappa=0}^{\mu-1}E_{n,l,\kappa}^{(1)} \\

\\

&  +h\displaystyle\sum_{l=0}^{n-1} \displaystyle\Bigg{(}\sum_{k=0}^{r-1}(\beta_{k}Q_{\mu-1}(z_{n}^{(\mu)},z_{l}^{(\mu-1)}-kh_{l}^{(\mu-1)})\delta_{h}(z_{l}^{(\mu-1)}-kh_{l}^{(\mu-1)})\\
& +\displaystyle\sum_{j=1}^{m}
\gamma_{j}Q_{\mu-1}(z_{n}^{(\mu)},z_{l}^{(\mu-1)}+c_{j}h_{l}^{(\mu-1)})\delta_{h}(z_{l}^{(\mu-1)}+c_{j}h_{l}^{(\mu-1)})\Bigg{)}
+h\displaystyle\sum_{l=0}^{n-1}E_{n,l,\mu-1}\\


& +h \displaystyle\sum_{\kappa=0}^{\mu-1} \displaystyle\sum_{l=0}^{N-1}\Bigg{(}\displaystyle\sum_{k=0}^{r-1}\beta_{k}Q_{\kappa}(z_{n}^{(\mu)},z_{l}^{(\kappa)}-kh_{l}^{(\kappa)})\delta_{h}(z_{l}^{(\kappa)}-kh_{l}^{(\kappa)})\\


& +\displaystyle\sum_{j=1}^{m}
\gamma_{j}Q_{\kappa}(z_{n}^{(\mu)},z_{l}^{(\kappa)}+c_{j}h_{l}^{(\kappa)})\delta_{h}(z_{l}^{(\kappa)}+c_{j}h_{l}^{(\kappa)})\Bigg{)}
+h\displaystyle\sum_{l=0}^{N-1}\displaystyle\sum_{\kappa=0}^{\mu-1}E_{n,l,v}^{(2)},
\end{array}
\end{equation}
Here, $E_{n,l,\mu}$, $E_{n,l,\mu-1}$, $E_{n,l,\kappa}^{(1)}$, and
$E_{n,l,\kappa}^{(2)}$ are the associated error terms. The assumption $c_{m}=1$, ensures that $z_{l-k}^{(\mu)}$, $z_{l-k}^{(\kappa)}$,
$z_{l,j}^{(\mu)}$, and $z_{l,j}^{(\kappa)}$ are collocation points for any choice of $k$ and $j$. As the defect function is zero at the collocation points, we get
$\delta(z_{l-k}^{(\mu)})=\delta(z_{l-k}^{(\kappa)})=\delta(z_{l,j}^{(\mu)})=\delta(z_{l,j}^{(\kappa)})=0$,
hence for $0\leq l <n \leq N-1$, $ 0\leq m \leq N-1$, we have
\begin{equation}\label{120}
 \begin{array}{lr}
\varepsilon(z_{n}^{(\mu)})=h\displaystyle\sum_{l=0}^{n-1}E_{n,l,\mu}+h\displaystyle\sum_{l=0}^{N-1}\displaystyle\sum_{\kappa=0}^{\mu-1}E_{n,l,\kappa}^{(1)}
+h\displaystyle\sum_{l=0}^{N-1}\displaystyle\sum_{\kappa=0}^{\mu-1}E_{n,l,\kappa}^{(2)}+h\displaystyle\sum_{l=0}^{n-1}E_{n,l,\mu-1}.
\end{array}
\end{equation}

 As a result, according to the
quadrature formula \eqref{118}, the order of
$\varepsilon(z_{n}^{(\mu)})$ is equal to the order of each
$E_{n,l,\mu}$, $E_{n,l,\mu-1}$, $E_{n,l,\kappa}^{(1)}$ and
$E_{n,l,\kappa}^{(2)}$. The condition \eqref{107} on the collocation
parameters  $\{c_{i}\}_{i=1}^{m}, $ ensures that the quadrature formula
\eqref{118} has the degree of $p-1$ and by employing the Peano
theorem with $p-1$, we can achieve

\begin{equation}
\max\limits_{1\leq n \leq N-1}\vert \varepsilon(z_{n}^{(\mu)})\vert
\leq Ch^{p},
\end{equation}
where the finite constant $C$ does not depend on $h$.
\end{proof}

%---------------------------------------------------
\section{ Linear Stability}

To analyze linear stability, we will assume the basic
test equation
\begin{equation}
  \label{1}
  y(z)=\left \{
  \begin{array}{lll}
    1+ \lambda\displaystyle\int_0^{z}y(\tau)d\tau+\lambda\displaystyle\int_0^{\Theta(z)}y(\tau)d\tau,& z\in[z_{0},Z],\\
    g (z),&z \in [\Theta(z_{0},z_{0}).
    \end{array}\right.
   \end{equation}
Now, we set
\begin{equation}
  \label{2}
  \begin{array}{llllll}
\displaystyle\int_0^{c_i}L_k(\tau)d\tau=\Omega_{ik},&\displaystyle\int_0^{c_i}\hat{L}(\tau)d\tau=\rho_{ij},\\
\displaystyle\int_0^{1}L_k(\tau)d\tau=\beta_{k},&\displaystyle\int_0^{1}\hat{L}(\tau)d\tau=\gamma_{j},\\
\displaystyle\int_{0}^1g(z_{l}^{(-1)}+\tau h_{l}^{(-1)})d\tau =\tilde{g}_l,&
\dst\int_{c_{i}}^{1}g(z_{l}^{(-1)}+\tau h_{l}^{(-1)})d\tau=\hat{g}_{n,i}
    \end{array}
\end{equation}
and we define %the vectors and matrices
\begin{equation}
  \label{3}
\begin{array}{lllll}
\mathbf{V}_n^{(\mu)}=\left[V_{n,1}^{(\mu)},...,V_{n,m}^{(\mu)}\right]^T,&\mathbf{y}_n^{(r,\mu)}=\left[y_n^{(\mu)},...,y_{n-r+1}^{(\mu)}\right]^T,\\
&\\
\mathbf{u}=\left[1,...,1\right]^T\in\rr^m,&\pmb\beta=\left[\beta_0,...,\beta_{r-1}\right]^T,\\
&\\
\pmb\gamma=\left[\gamma_1,...,\gamma_m\right]^T,&
\pmb{\hat{L}}(1)=\left[\hat{L}_1(1),...,\hat{L}_m(1)\right]^T,\\
&\\
\pmb{L}(1)=\left[L_0(1),...,L_{r-1}(1)\right]^T,&\pmb{L}(0)=\left[L_0(0),...,L_{r-1}(0)\right]^T,\\
&\\
\pmb{\hat{\Omega}}=\left(\Omega_{ik}\right)\in\rr^{m\times r},
&\pmb{\rho}=\left(\rho_{ij}\right)\in\rr^{m\times m},\\
 &\\
 \pmb{\mathcal{A}}(z)=\left[\begin{array}{c}
                   -\pmb{\hat{L}}^T(z) \\\hline
                                         \mathbf{0}_{r\times m}
                                                      \end{array}
        \right],~~~z=0,1,&\pmb{\mathcal{C}}(z)=\left[\begin{array}{c|c}
                                 1 & -\pmb{L}^T(z) \\\hline
                                 \mathbf{0}_{r\times1}&\mathbf{I}_r
                               \end{array}
                  \right],~~~z=0,1,
                  &\\
                  &\\
                \pmb{\mathcal{B}}=\left[\begin{array}{c|c}
                                 \mathbf{0}_{1\times r} & 0 \\\hline
                                 \mathbf{I}_r & \mathbf{0}_{r\times 1}
                                                      \end{array}
        \right],&\tilde{\pmb G}_{n}=\displaystyle\int_{0}^1 g(t_{n}^{(-1)}+\tau h_{n}^{(-1)})d\tau ,\\
        \end{array}
\end{equation}
and
$$
\pmb{\hat{G}}_{n}
=\left[\displaystyle\int_{c_1}^1g(z_{n}^{(-1)}+\tau h_{n}^{(-1)})d\tau,
 \ldots,\displaystyle\int_{c_m}^1\phi(z_{n}^{(-1)}+\tau h_{n}^{(-1)})\right]^T.
$$

\begin{theorem}
\label{thst1} For the test equation \eqref{1},  by using the exact
MSCM \cite{darania-sutudeh}, we  have the
following recurrence relation

\begin{equation} \label{5}
 \pmb{\mathcal{Y}}_{n} =\pmb{\mathcal{T}}(\zeta)\pmb{\mathcal{Y}}_{n-1}+\pmb{\mathcal{N}}(\zeta),
     \end{equation}
where
\begin{equation}\label{6}
 \pmb{\mathcal{Y}}=\left[
 \begin{array}{lllll}
y_{n+1}^{(\mu)}\\
\mathbf{y}_n^{(r,\mu)}\\
\mathbf{V}_{n}^{(\mu)}\\
y_{n}^{(\mu-1)}\\
\mathbf{y}_{n}^{(r,\mu-1)}\\
\mathbf{V}_{n}^{(\mu-1)}\\
  \end{array}\right],
   \end{equation}
and the stability matrix is $\pmb{\mathcal{T}}(\zeta)\in \rr^{(p+r+3)^2}$,  which is provided by
\begin{equation}
  \label{7}
\pmb{\mathcal{T}}(\zeta)=\left[\pmb{\mathcal{P}}(\zeta)\right]^{-1}\pmb{\mathcal{M}}(\zeta),~~\pmb{\mathcal{N}}(\zeta)=\left[\pmb{\mathcal{P}}(\zeta)\right]^{-1}\pmb{\mathcal{J}}(\zeta),
     \end{equation}
with $\zeta=\lambda h$ and
\begin{equation}
  \label{8}
  \pmb{\mathcal{P}}(\zeta)=\left[
  \begin{array}{c|c|c|cc|c}
 \begin{array}{c|c}\mathbf{0} & -\zeta\pmb{\hat{\Omega}}\end{array}&\mathbf{I}_m-\zeta\pmb\rho&\begin{array}{c|c}\mathbf{0}&-\zeta\pmb{\hat{\Omega}}\end{array}&-\zeta\pmb\rho\\\hline
 \begin{array}{c|c} \mathbf{0} & -\zeta\pmb{\hat{\Omega}}\end{array}&\mathbf{I}_m-\zeta \pmb\rho&\begin{array}{r|r}\mathbf{0}&~~~~\mathbf{0}\end{array}&\mathbf{0}\\\hline
 \pmb{\mathcal{C}}(1)&\pmb{\mathcal{A}}(1)&\mathbf{0}&\mathbf{0}\\\hline
\mathbf{0}&\mathbf{0}&\pmb{\mathcal{C}}(0)&\pmb{\mathcal{A}}(0)
  \end{array} \right],
     \end{equation}
\begin{equation}
  \label{9}
  \pmb{\mathcal{M}}(\zeta)=\left[
   \begin{array}{c|c|c|cc|c}
 \begin{array}{c|c} \mathbf{0} &\zeta(\pmb u\pmb\beta^T-\pmb{\hat{\Omega}})\end{array}&\mathbf{I}_m+\zeta(\pmb u\pmb\gamma^T-\pmb\rho)&\begin{array}{c|c}\mathbf{0}&\zeta(\pmb u\pmb\beta^T-\pmb{{\hat{\Omega}})}\end{array}&\zeta(\pmb u\pmb\gamma^T-\pmb\rho)
 \\\hline
  \begin{array}{c|c} \mathbf{0} & \zeta(\pmb u\pmb\beta^T-\pmb{\hat{\Omega}})\end{array}&\mathbf{I}_m+\zeta(\pmb u\pmb\gamma^T-\pmb\rho)&\begin{array}{c|c}\mathbf{0}&~~~~~~~~\mathbf{0}~~~~~\end{array}&\mathbf{0}\\\hline
 \pmb{\mathcal{B}}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\hline
\mathbf{0}&\mathbf{0}&\pmb{\mathcal{B}}&\mathbf{0}
  \end{array} \right],
     \end{equation}

\begin{equation}
  \label{10}
\pmb{\mathcal{J}}(\zeta)=-\zeta\left[
\begin{array}{c}
 \mathbf{0} \\
\hat{\pmb G}_{n-1}-\hat{\pmb G}_{n}+\pmb u\tilde{\pmb G}_{n}\\
 \mathbf{0}\\
 \mathbf{0} \\
 \end{array} \right].
\end{equation}
%and $ X_{n-1,n}=\hat{\pmb G}_{n-1}-\hat{\pmb G}_{n}$.

     \end{theorem}
\begin{proof}
The collocation equation of \eqref{1} is
\begin{align}\label{11}
v_{h}(z)=1+\lambda\displaystyle\int_{z_{0}}^{z}v_{h}(\tau )d\tau+\lambda\displaystyle\int_{z_{0}}^{\Theta(z)}v_{h}(\tau)d\tau,
\end{align}

by substituting $z=z_{n,i}^{(\mu)}$, in the equation \eqref{11}, we
obtain
\begin{equation}\label{12}
\begin{array}{rlll}


v_{h}(z_{n,i}^{(\mu)})=&1+\lambda\displaystyle\int_{z_{0}}^{z_{n,i}^{(\mu)}}v_{h}(\tau )d\tau
+\lambda\displaystyle\int_{z_{0}}^{\Theta(z_{n,i}^{(\mu)})}v_{h}(\tau)d\tau\\
=&1+\lambda(\mathcal{V}v_{h})(z_{n,i}^{(\mu)})+\lambda(\mathcal{V}_{\Theta}v_{h})(z_{n,i}^{(\mu)}),
\end{array}
\end{equation}

where
\begin{equation}\label{13}
\begin{array}{lr}
(\mathcal{V}v_{h})(z_{n,i}^{(\mu)})=\displaystyle\sum_{\eta=0}^{\mu-1}\displaystyle\sum_{l=0}^{N-1}h_{l}^{(\eta)}
\displaystyle\int_{0}^{1}v_{h}(z_{l}^{(\eta)}+\tau h_{l}^{(\eta)})d\tau \\
\qquad \qquad \qquad
+\displaystyle\sum_{l=0}^{n-1}h_{l}^{(\mu)}\int_{0}^{1}v_{h}(z_{l}^{(\mu)}+\tau h_{l}^{(\mu)})d\tau +h_{n}^{(\mu)}\int_{0}^{c_{i}}v_{h}(z_{n}^{(\mu)}+\tau h_{n}^{(\mu)})d\tau,

\end{array}
\end{equation}
and
\begin{equation}\label{14}
(\mathcal{V_{\theta}}v_{h})(z_{n,i}^{(\mu)})=\left\{
\begin{array}{lr}
 \displaystyle{\sum_{\eta=0}^{\mu-2}}\sum_{l=0}^{N-1}h_{l}^{(\eta)}\int_{0}^{1}v_{h}(z_{l}^{(\eta)}+\tau h_{l}^{(\eta)})d\tau \\
+\displaystyle{\sum_{l=0}^{n-1}h_{l}^{(\mu-1)}\int_{0}^{1}}v_{h}(z_{l}^{(\mu-1)}
+\tau h_{l}^{(\mu-1)})d\tau \\
+h_{n}^{(\mu-1)}\displaystyle\int_{0}^{c_{i}}v_{h}(z_{n}^{(\mu-1)}+\tau h_{n}^{(\mu-1)})d\tau ,
\quad \Theta(z_{n,i}^{(\mu)}>z_{0},~\mu=1,\ldots
,M,\\\\
 -h_{n}^{(-1)}\displaystyle\int_{c_{i}}^{1}g(z_{n}^{(-1)}
 +\tau h_{n}^{(-1)}))d\tau \\
-\displaystyle\sum_{l=n+1}^{N-1}h_{l}^{(-1)}\displaystyle\int_{0}^{1}g(t_{l}^{(-1)}+\tau h_{l}^{(-1)})d\tau ,
 \qquad \Theta(z_{n,i}^{(0)}) \leq z_{0}.
\end{array}\right.
\end{equation}
In equations \eqref{13} and \eqref{14}, we assume that
$h=\max_{\kappa,l}\limits h_l^{(\kappa)},~\kappa=0, \ldots, \mu$, and $\zeta=h\lambda$. Then by
using the equation \eqref{14*}, we have
\begin{equation}\label{19}
\pmb V_{n}^{(\mu)}=\left\{
\begin{array}{lr}
\pmb u+\zeta\bigg\{\pmb
u\bigg[\displaystyle\sum_{\eta=0}^{\mu-1}\displaystyle\sum_{l=0}^{N-1}(\pmb\beta^{T}\pmb
y_{l}^{(r,\eta)}+\pmb\gamma^{T}V_{l}^{(\eta)})\\+
\displaystyle\sum_{l=0}^{n-1}(\pmb\beta^{T}\pmb y_{l}^{(r,\mu)}
 +\pmb\gamma^{T}V_{l}^{(\mu)})
  +\displaystyle\sum_{\eta=0}^{\mu-2}\displaystyle\sum_{l=0}^{N-1}(\pmb\beta^{T}\pmb y_{l}^{(r,\eta)}
 +\pmb\gamma^{T}\pmb V_{l}^{(\eta)})\\
 +\displaystyle\sum_{l=0}^{n-1}(\pmb\beta^{T}\pmb y_{l}^{(r,\mu-1)}+\pmb\gamma^{T}\pmb V_{l}^{(\mu-1)})\bigg]\\
 +\pmb{\hat{\Omega}} (\pmb y_{n}^{(r,\mu)}+ \pmb y_{n}^{(r,\mu-1)})\\+\pmb\rho(\pmb V_{n}^{(\mu)}+\pmb V_{n}^{(\mu-1)})
 \bigg\}, \quad  \quad  \quad \Theta(z_{n,i}^{(\mu)})>z_{0},~\mu=1,\ldots,M,\\\\

\pmb u+\zeta\bigg {\{ }\pmb
u\bigg[\displaystyle\sum_{\eta=0}^{\mu-1}\displaystyle\sum_{l=0}^{N-1}(\pmb\beta^{T}\pmb
y_{l}^{(r,\eta)})+\pmb\gamma^{T}V_{l}^{(\eta)} )\\+
\displaystyle\sum_{l=0}^{n-1}(\pmb\beta^{T}\pmb y_{l}^{(r,\mu)}
 +\pmb\gamma^{T}V_{l}^{(\mu)})- \displaystyle\sum_{l=n+1}^{N-1}\tilde{G}_{l}\bigg]
 +\pmb {\hat{\Omega}} \pmb y_{n}^{(r,\mu)}+\pmb\rho\pmb V_{n}^{(\mu)}\\
-\pmb{\hat{G}}_{n}\bigg{ \} },  \quad  \quad \quad  \Theta(z_{n,i}^{(0)})\leq z_{0},
 \end{array}\right.
\end{equation}
The computation of the difference
$\mathbf{V}^{(\mu)}_n-\mathbf{V}^{(\mu)}_{n-1}$ by inserting the
\eqref{19} for both expressions $\mathbf{V}^{(\mu)}_n$ and
$\mathbf{V}^{(\mu)}_{n-1}$, when $\Theta(t_{n,i}^{(\mu)})
>z_{0}$ and $\Theta(t_{n,i}^{(\mu)}) \leq z_{0}$, leads to

\begin{equation}\label{31}
\pmb Q_{1}(\zeta) \begin{bmatrix}
\pmb V_{n}^{(\mu)}\\
\pmb V_{n}^{(\mu-1)}\\
\pmb y_{n}^{(r,\mu)}\\
\pmb y_{n}^{(r,\mu-1)}
\end{bmatrix}
= \pmb M_{1}(\zeta)
\begin{bmatrix}
\pmb V_{n-1}^{(\mu)}\\
\pmb V_{n-1}^{(r,\mu-1)}\\
\pmb y_{n-1}^{(r,\mu)}\\
\pmb y_{n-1}^{(r,\mu-1)}
\end{bmatrix}
+ \pmb N_1(\zeta),
\end{equation}
where
\begin{equation}\label{28-1}
\pmb Q_{1}(\zeta)=
\begin{bmatrix}
I_{m}-\zeta\pmb \rho & -\zeta\pmb \rho & -\zeta\pmb \rho & -\zeta\pmb\Omega\\

I_{m}-\zeta\pmb \rho & 0 & -\zeta\pmb{\hat{\Omega}} & 0,
\end{bmatrix},
\end{equation}

\begin{equation}\label{29}
\pmb M_{1}(\zeta)=
\begin{bmatrix}
I_{m}+\zeta(\pmb\gamma^{T}-\pmb\rho) &\zeta(\pmb\gamma^{T}-\pmb\rho) & \zeta(\pmb\beta^{T}-\pmb{\hat{\Omega}}) & \zeta(\pmb\beta^{T}-\pmb{\hat{\Omega}})\\
I_{m}+\zeta(\pmb\gamma^{T}-\pmb\rho) & 0 & \zeta(\pmb \beta^{T}-\pmb{\hat{\Omega}})
& 0,
\end{bmatrix},
\end{equation}

\begin{equation}\label{30}
\pmb N_1(\zeta)= -\zeta\pmb u
\begin{bmatrix}
\pmb 0\\
\hat{\pmb G}_{n-1}-\hat{\pmb G}_{n}
\end{bmatrix}.
\end{equation}
Using the notations in this section and based on the provided equation
\[
 y(z^{(\mu)}_n+\tau h)=\displaystyle\sum_{k=0}^{r-1}L_k(\tau )y^{(\mu)}_{n-k}+\displaystyle\sum_{j=1}^{m}\hat{L}(\tau )Y^{(\mu)}_{n,j},~~~\tau \in[0,1],
\]
we have
\begin{equation}
 \label{es9}
 y^{(\mu)}_{n+1}=\pmb{L}^T(1)\mathbf{y}^{(r,\mu)}_n+\pmb{\hat{L}}^T(1)\mathbf{V}^{(\mu)}_n,
     \end{equation}
and
\begin{equation}
 \label{es9-1}
 y^{(\mu-1)}_{n}=\pmb{L}^T(0)\mathbf{y}^{(r,\mu-1)}_n+\pmb{\hat{L}}^T(0)\mathbf{V}^{(\mu-1)}_n,
     \end{equation}
We can express it in the following matrix forms:
%We can express it in the following matrices  forms:
\begin{equation}
  \label{es10}
\pmb{\mathcal{C}}(1)\left[
  \begin{array}{lllll}
y^{(\mu)}_{n+1}\\
\mathbf{y}_n^{(r,\mu)}
   \end{array} \right]+\pmb{\mathcal{A}}(1)\mathbf{V}^{(\mu)}_n=\pmb{\mathcal{B}}\left[
  \begin{array}{lllll}
y^{(\mu)}_{n}\\
\mathbf{y}_{n-1}^{(r,\mu)}
   \end{array} \right],
     \end{equation}
\begin{equation}
  \label{es10-1}
\pmb{\mathcal{C}}(0)\left[
  \begin{array}{lllll}
y^{(\mu-1)}_{n}\\
\mathbf{y}^{(r,\mu-1)}_{n}
   \end{array} \right]+\pmb{\mathcal{A}}(0)\mathbf{V}^{(\mu-1)}_{n}=\pmb{\mathcal{B}}\left[
  \begin{array}{lllll}
y^{(\mu-1)}_{n-1}\\
\mathbf{y}^{(r,\mu-1)}_{n-1}
   \end{array} \right],
     \end{equation}
where $\pmb{\mathcal{A}}(z),\pmb{\mathcal{C}}(z),~z=0,1$  and $\pmb{\mathcal{B}}$ are
given by \eqref{3}.



 Now, we conclude from
\eqref{es10}, \eqref{es10-1}  and \eqref{31}-\eqref{30} that
\begin{equation}\label{36}
\begin{array}{lr}
\pmb{\mathcal{ P}}(\zeta)
\begin{bmatrix}
\pmb y_{n+1}^{(\mu)}\\
\pmb y_{n}^{(r,\mu)}\\
\pmb V_{n}^{(\mu)}\\
\pmb y_{n}^{(\mu-1)}\\
\pmb y_{n}^{(r,\mu-1)}\\
\pmb V_{n}^{(\mu-1)}
\end{bmatrix}
= \pmb{\mathcal{ M}}(\zeta)
\begin{bmatrix}
\pmb y_{n}^{(\mu)}\\
\pmb y_{n-1}^{(r,\mu)}\\
\pmb V_{n-1}^{(\mu)}\\
\pmb y_{n-1}^{(\mu-1)}\\
\pmb y_{n-1}^{(r,\mu-1)}\\
\pmb V_{n-1}^{(\mu-1)}
\end{bmatrix}
+ \pmb J(\zeta)
\end{array}
\end{equation}
where
\begin{equation*}
 \pmb{\mathcal{P}}(\zeta)=\left[
  \begin{array}{c|c|c|cc|c}
 \begin{array}{c|c}\mathbf{0} & -\zeta\pmb\Omega\end{array}&\mathbf{I}_m-\zeta\pmb\rho&\begin{array}{c|c}\mathbf{0}&-\zeta\pmb\Omega\end{array}&-\zeta\pmb\rho\\\hline
 \begin{array}{c|c} \mathbf{0} & -\zeta\pmb\Omega\end{array}&\mathbf{I}_m-\zeta\pmb\rho&\begin{array}{r|r}\mathbf{0}&~~~~\mathbf{0}\end{array}&\mathbf{0}\\\hline
\pmb{\mathcal{C}}(1)&\pmb{\mathcal{A}}(1)&\mathbf{0}&\mathbf{0}\\\hline
\mathbf{0}&\mathbf{0}&\pmb{\mathcal{C}}(0)&\pmb{\mathcal{A}}(0)
  \end{array} \right],
     \end{equation*}

\begin{equation*}
\pmb{\mathcal{M}}(\zeta)=\left[
   \begin{array}{c|c|c|cc|c}
 \begin{array}{c|c} \mathbf{0} & \zeta(\pmb u\pmb\beta^T-\pmb\Omega)\end{array}&\mathbf{I}_m+\zeta(\pmb u\pmb\gamma^T-\pmb\rho)&\begin{array}{c|c}\mathbf{0}&\zeta(\pmb u\pmb\beta^T-\pmb\Omega)\end{array}&\zeta(\pmb u\pmb\gamma^T-\pmb\rho)
 \\\hline
  \begin{array}{c|c} \mathbf{0} &\zeta(\pmb u\pmb\beta^T-\pmb\Omega)\end{array}&\mathbf{I}_m+\zeta(\pmb u\pmb\gamma^T-\pmb\rho)&\begin{array}{c|c}\mathbf{0}&~~~~~~~~\mathbf{0}~~~~~\end{array}&\mathbf{0}\\\hline
\pmb{\mathcal{B}}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\hline
\mathbf{0}&\mathbf{0}&\pmb{\mathcal{B}}&\mathbf{0}
  \end{array} \right],
     \end{equation*}

\begin{equation*}
\pmb{\mathcal{J}}(\zeta)=-\zeta\left[
\begin{array}{c}
 \mathbf{0} \\
\hat{\pmb G}_{n-1}-\hat{\pmb G}_{n}+\pmb u\tilde{\pmb G}_{n}\\
 \mathbf{0}\\
 \mathbf{0} \\
 \end{array} \right],
\end{equation*}
or more concisely, as
\begin{equation}\label{38}
 \mathbf{\mathcal{Y}}_{n}=\pmb{\mathcal{T}}(\zeta) \mathbf{\mathcal{Y}}_{n-1} +\pmb{\mathcal{N}}(\zeta),
\end{equation}
where
 $\pmb{\mathcal{ T}}(\zeta)=\pmb{\mathcal{ P}}^{-1}(\zeta) \pmb{\mathcal{ M}}(\zeta),~\pmb{\mathcal{N}}(\zeta)=\pmb{\mathcal{ P}}^{-1}(\zeta)
\pmb J(\zeta)$ are known, thus concluding the proof.
%are known and thus concluding the proof.
\end{proof}

\begin{theorem}
\label{thst2} Let the assumption of Theorem \ref{thst1}
holds and $\|\pmb{\mathcal{T}}(\zeta)\|\leq1$. Consequently, the equation \eqref{38} is stable.
\end{theorem}

\begin{proof}
The proof is derived by making a simple adjustment to Theorem 3.3, as presented in \cite{darania-pishbin}.
\end{proof}

\begin{theorem}
\label{thst3} In the test equation \eqref{1}, by using the
discretized  MSCM \cite{darania-sutudeh}, we
get the following recurrence equation
\[
\pmb{\mathcal{Y}}_{n} =\pmb{\mathcal{{\tilde{T}}}}(\zeta)\pmb{\mathcal{Y}}_{n-1}+\pmb{\mathcal{ {\tilde{N}}}}(\zeta),
     \]
where  $\pmb{\mathcal{A}}(.),~\pmb{\mathcal{C}}(.), ~\pmb{\mathcal{J}}(.)$, and
$\pmb{\mathcal{B}}$  are already defined.
 Also,  the matrix of stability  $\pmb{\mathcal{\tilde{T}}}(\zeta)$, is provided by
$$
  \pmb{\mathcal{\tilde{T}}}(\zeta)=\left[\pmb{\mathcal{\tilde{P}}}(z)\right]^{-1}\pmb{\mathcal{\tilde{M}}}(\zeta)
  ,~~\pmb{\mathcal{\tilde{N}}}(\zeta)=\left[\pmb{\mathcal{\tilde{P}}}(\zeta)\right]^{-1}\pmb{\mathcal{J}}(\zeta),
    $$
and
$$
   \pmb{\mathcal{\tilde{P}}}(\zeta)=\left[
  \begin{array}{c|c|c|cc|c}
 \begin{array}{c|c}\mathbf{0} & -\zeta\tilde{\pmb\Omega}\end{array}&\mathbf{I}_m-\zeta\tilde{\pmb\rho}&\begin{array}{c|c}\mathbf{0}&-\zeta\tilde{\pmb\Omega}\end{array}&-\zeta\tilde{\pmb\rho}\\\hline
 \begin{array}{c|c} \mathbf{0} & -\zeta\tilde{\pmb\Omega}\end{array}&\mathbf{I}_m-\zeta\tilde{\pmb\rho}&\begin{array}{c|c}\mathbf{0}&~~~~\mathbf{0}\end{array}&\mathbf{0}\\\hline
 \pmb{\mathcal{C}}(1)&\pmb{\mathcal{A}}(1)&\mathbf{0}&\mathbf{0}\\\hline
\mathbf{0}&\mathbf{0}&\pmb{\mathcal{C}}(0)&\pmb{\mathcal{A}}(0)
  \end{array} \right],
    $$

$$
    \pmb{\mathcal{\tilde{M}}}(\zeta)=\left[
   \begin{array}{c|c|c|cc|c}
 \begin{array}{c|c} \mathbf{0} & \zeta(\pmb u\tilde{\pmb\beta}^T-\tilde{\pmb\Omega})\end{array}&\mathbf{I}_m+\zeta(\pmb u\tilde{\pmb\gamma}^T-\tilde{\pmb\rho})&\begin{array}{c|c}
 \mathbf{0}&\zeta(\pmb u\tilde{\pmb\beta}^T-\tilde{\pmb\Omega})\end{array}&\zeta(\pmb u\tilde{\pmb\gamma}^T-\tilde{\pmb\rho})\\\hline
  \begin{array}{c|c} \mathbf{0} & \zeta(\pmb u\tilde{\pmb\beta}^T-\tilde{\pmb\Omega})\end{array}&\mathbf{I}_m+\zeta(\pmb u\tilde{\pmb\gamma}^T
  -\tilde{\pmb\rho})&\begin{array}{l|l}\mathbf{0}&~~~~\mathbf{0}~~~~~~~~~\end{array}&\mathbf{0}\\\hline
\pmb{\mathcal{B}}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\hline
\mathbf{0}&\mathbf{0}&\pmb{\mathcal{B}}&\mathbf{0}
  \end{array} \right],
   $$
and

$$
  \begin{array}{lllll}
   \tilde{\Omega}_{jk}=\displaystyle\sum_{l=1}^{\mu_0}\limits \omega_{jl}L_{k}(d_{jl}), &
   \tilde{\rho}_{ik}=\displaystyle\sum_{l=1}^{\mu_0}\limits\omega_{il}\hat{L}_{k}(d_{il}),&
   \tilde{\beta}_{k}=\displaystyle\sum_{l=1}^{\mu_1}\limits \omega_lL_{k}(d_l),&
   \tilde{\gamma}_{j}=\displaystyle\sum_{l=1}^{\mu_1}\limits\omega_l\hat{L}_{j}(d_l),\\
   &&&\\
   \tilde{\pmb\beta}=[\tilde{\beta}_0,...,\tilde{\beta}_{r-1}]^T,&
   \tilde{\pmb\gamma}=[\tilde{\gamma}_1,...,\tilde{\gamma}_m]^T,&
   \tilde{\pmb\Omega}=(\tilde{\Omega}_{ik})\in\rr^{m\times r},&
   \tilde{\pmb\rho}=(\tilde{\rho}_{ij})\in\rr^{m\times m}.
\end{array}
$$
\end{theorem}
\begin{proof} The proof follows immediately from the consequence of Theorem \ref{thst1}.
%It resembles the proof of Theorem \ref{thst1}.
\end{proof}


The stability function of the method with regard to \eqref{1} is given by
\begin{equation}
 \label{st1}
Q(x,\zeta)=\det (x\mathbf{I}_{p+r+3}-\pmb{\mathcal{T}}(\zeta)).
 \end{equation}
Additionally, the stability properties of the exact MSCM will be analyzed with the polynomial acquired by multiplying the stability
 function \eqref{st1} by its denominator. The resulting polynomial will be denoted by the same
  $ Q(x,\zeta)$. Moreover, the area of absolute stability of the methods is defined by
$$
\pmb{\mathcal{S}}:=\{\zeta\in\Bbb{C}: |x_l(\zeta)|<1,~~l=1,2,...,p+r+3\},
$$
which is represented by the roots of the polynomial $p(x,\zeta)$, namely $\{x_l\}_{l=1}^{p+r+3}$.
For this method, we have
\begin{equation}
 \label{st2}
Q(x,\zeta)=\sum_{l=0}^{p+r+3}Q_l(\zeta)x^l,
 \end{equation}
 where $Q_l(\zeta ),~ l=0,1,...,p+r+3$, refers to polynomials of degree $2m$ or lower.
 We will use the boundary locus method \cite{Liniger} to achieve the area of absolute stability by substituting $ x = e^{i\theta}$, the roots of \eqref{st2} define the area of stability.
%****************************************************************\\

\subsection{Stability area}
The present section will demonstrate the analytical results achieved in Section 3, using the following instances.
Also, the stability polynomial for MSCM can be represented by the form \eqref{st2}. It should be noted that Fig. 1,  illustrates the stability area for the MSCM with $c = (0.31,0.75,1)$ and $r =m = 3$.  Meanwhile, Figures 2, and 3  depict the stability area for the MSCM with  $ r = 4, m =
2, c =( 0.7,  1)$ and $m = 2, r = 3$ including the super  convergence collocation parameters $c = (\frac{21}{38},1)$, respectively.

\begin{rem}

\label{rems1}
The order of applied quadrature rules in the discretized MSCM is shown to be at least the same order for the MSCM  in Section 2. Since the polynomials $L_k(s), k = 0, 1, ..., r-1$ and $\hat{L}(s), j = 1,
2,..., m,$ possess the degree of $2m+ 2r+ 2$, the quadrature rules will be exact for them. Moreover, for the MSCM, we have $\pmb{\mathcal{T}} (\zeta) =\pmb{\mathcal{ \tilde{T}}}(\zeta)$. Stability areas are displayed in Figures 1, 2, and 3, and will remain unchanged for the discretized cases.
\end{rem}
\begin{figure}[h!]\hspace{3cm}
\label{figer1}
\begin{center}
\scalebox{.5} {\includegraphics{Fig1.eps}}\\
\caption{{\small
The region of stability for the MSCM (convergence).
}}
\end{center}
\end{figure}
\begin{figure}[h!]\hspace{3cm}
\label{figer2}
\begin{center}
\scalebox{.5} {\includegraphics{Fig2.eps}}\\
\caption{{\small The region of stability for the MSCM  (convergence).}}
\end{center}
\end{figure}
\begin{figure}[h!]\hspace{3cm}
\label{figer3}
\begin{center}
\scalebox{.5} {\includegraphics{Fig3.eps}}\\
\caption{{\small
The region of stability for the MSCM   (super-convergence).}}
\end{center}
\end{figure}


\goodbreak

\section{ Numerical examples }
 In this segment, some examples will be evaluated to validate our theoretical results. We are considering the MSCM  with $ r=3$ and $ m=2$, as shown in Examples \ref{ex1} and \ref{ex2}, which were selected from the collocation abscissas of the system \eqref{107}, with  $(c_1, c_2)=(\frac{21}{38},1)$. The method has an order of $p = 6$ (super-convergence).

\begin{exam}\cite{darania-sutudeh}
\label{ex1}
 Let us consider the non-linear delay integral equations as follows:
$$y(z)= \left \{
  \begin{array}{lll}
    \phi(z)+\displaystyle\int_{\frac{1}{4}}^{z}(\tau ^2+z+1)y^2(\tau)d\tau+\int_{\frac{1}{4}}^{\frac{1}{2}z}(2\tau+z^2+4)y^2(\tau)d\tau,& z\in (\frac{1}{4},1],\\
    &\\
    z+ \sin z ,&z\in[\frac{1}{8}, \frac{1}{4}].
    \end{array}\right.$$
here,
 $\phi (z)$ is such that  the  exact solution is
$y(z)=z+ \sin z.$
\end{exam}

In Example \ref{ex1}, we  assume that   $\xi_{2}=Z=1$ , then
$$\xi_{\mu}=(\frac{1}{2})^{2-\mu},~~~ \mu=0,1.$$

The highest (maximum) errors at the grid points for $r=3$ and $m = 2$
versus $N$ are presented in Table 1. The rate of convergence that was observed is computed from the maximum errors at the grid points, which corroborate the results of the convergence in theorem \ref{th4}.

Observe  that, the convergence order for $c=(0.7,1)$, is $p_{c}=m+r$, and the super-convergence order for $ c=(\frac{21}{38},1)$, is $p_{sc}=2m+r-1$.

\begin{center}
\begin{tabular}{c c c c c c c c c|}
\hline & &  $  c=(\frac{21}{38},1)$  & & &   $  c=(0.7,1)$   & &   \\
 \cline{3-3}    \cline{6-6} {\small{\it N}} & &
{\small $||y-v_h||_{\infty}$}~~~~~~$p_{sc}$ &&& {\small $||y-v_h||_{\infty}$}~~~~~~$p_{c}$   &   \\
\hline \vspace{-.22cm}\\
 4& &$1.62\times10^{-8}$~~~~~4.95 && &$5.461\times10^{-8}$~~~~~2.72 \\
 8& &$1.07\times10^{-9}$ ~~~~~5.46&&& $8.28\times10^{-9}$~~~~~ 4.25\\
16& &$2.43\times10^{-11}$~~~~~5.73 &&&$ 4.33   \times10^{-10}$~~~~~4.68\\
32& &$4.57\times10^{-13}$~~~~~~~~-&&&$  1.68 \times10^{-11}$ ~~~~~-\\
\hline\\
\end{tabular}~~~~~~~~~~~~~~~~~
\\ \scriptsize{\textbf{Table 1}.
The highest errors $\|y-v_h\|_{\infty}$  resulting from $m=2$ and $r=3$
 in Example \ref{ex1}.
}
\end{center}


\begin{exam}\cite{darania-sutudeh}
\label{ex2} Consider the linear delay integral equations with
the non-linear lag function $\Theta(t)$ as:
$$y(z)= \left \{
  \begin{array}{lll}
    z+\displaystyle\int_{\frac{1}{4}}^{\frac{1}{2}z^2}\tau y(\tau)d\tau,&z\in (\frac{1}{4},1],\\
    &\\
    e^z ,&z\in[\frac{1}{32}, \frac{1}{4}].
    \end{array}\right.$$
where
  the  exact solution is
$$y(z)= \left \{
  \begin{array}{lll}
   e^z,&z\in[\frac{1}{32}, \frac{1}{4}],\\
  \frac{3}{4} e^{\frac{1}{4}} + z +
 \frac{1}{2}e^{\frac{z^2}{2}} (-2 + z^2), & z\in(\frac{1}{4}, \frac{1}{\sqrt{2}}],\\
z +\frac{1}{384} (-2 + 756 e^{\frac{1}{32}} + 16z^6 +
    48 e^{\frac{z^4}{8}} (-16 +zt^4) +
    9 e^{\frac{1}{4}} (-1 + 4 z^4)) ,&z\in(\frac{1}{\sqrt{2}}, 2^{\frac{1}{4}}].
    \end{array}\right.$$
\end{exam}
In Example \ref{ex2}, we identify that the solution of the delay equation is discontinuous at $z_0=\frac{1}{4}$. Hence, the primary discontinuity points $\xi_{\mu}$ can be attributed to the formula as follows:

$$\Theta(\xi_{\mu})=\frac{1}{2}(\xi_{\mu})^2 =\xi_{\mu-1}, ~~~\mu=1,\cdots~~~;\xi_{0}=z_0=\frac{1}{4},  $$
and $Z=1\in(\xi_{1}, \xi_{2}].$ Also,  the images $\{ \tilde{c}_{i}
\}$ of the collocation parameters  $\{c_{i} \}$ can be obtained from
the following relation
$$\Theta(z_{n}^{(\mu)}+c_{i}h_{n}^{(\mu)})=\frac{1}{2}(z_{n}^{(\mu)}+c_{i}h_{n}^{(\mu)})^2=z_{n}^{(\mu-1)}+\tilde{c}_{i}h_{n}^{(\mu-1)},\qquad i=1, \ldots , m.$$
% The maximum errors have been displayed for different values of $N$ in Table 2.
Table 2 shows the highest errors for some values of $N$.  These results are consistent with the theoretical prediction.
\begin{center}
\begin{tabular}{c c c c c c c c c|}
\hline & &  $  c=(\frac{21}{38},1)$  & & &   $  c=(0.7,1)$   & &   \\
 \cline{3-3}    \cline{6-6} {\small{\it N}} & &
{\small $||y-v_h||_{\infty}$}~~~~~~$p_{sc}$  &&& {\small $||y-v_h||_{\infty}$}~~~~~~$p_{c}$    &   \\
\hline \vspace{-.22cm}\\
 4& &$2.62\times10^{-8}$~~~~~4.60 &&&$3.16 \times10^{-8} $~~~~~3.99  \\
 8& &$1.96\times10^{-9}$ ~~~~~5.42& &&$1.99\times10^{-9} $~~~~~ 4.75 \\
16& &$4.55\times10^{-12}$~~~~~5.65 & &&$7.38\times10^{-11} $~~~~  4.90\\
32& &$9.04\times10^{-14}$~~~~~~~~-& &&$ $~~~~~-  \\
\hline\\
\end{tabular}~~~~~~~~~~~~~~~~~
\\ \scriptsize{\textbf{Table 2}.The highest errors $\|y-v_h\|_{\infty}$ resulting from $m=2$ and $r=3$
 in Example \ref{ex1}.
}
\end{center}\goodbreak


%
\section{Conclusion}
This paper describes a multistep collocation technique for approximating the solutions of delay integral equations. In addition, we present a general analysis of the super-convergence and linear stability of this method. Furthermore, the multistep collocation methods exhibit a uniform order of $m+r$ for all choices of collocation parameters \cite{darania-sutudeh} and a local super-convergence order of $2m+r-1$ in the mesh points of collocation parameters. In order to evaluate the efficiency of the preferred method, various instances are conducted.


\goodbreak
%\section* {Acknowledgments}
%
%The authors are very grateful to the Editor and anonymous reviewers
%for their valuable comments and suggestions.

% ----------------  ------------------------------------------------------
\goodbreak


\goodbreak
 \begin{thebibliography}{99}
%% ----------------------------------------------------------------------
\bibitem{darania-sutudeh} P. Darania and F. Sutudeh Maram, Numerical analysis of a high order method for nonlinear delay
integral equations, Journal of Computational and Applied Mathematics
374 (2020) 112738.
% ----------------------------------------------------------------------

\bibitem{darania-tabriz} P. Darania, super-convergence analysis of multi-step collocation method for delay
 Volterra integral equations, Computational Methods for
Differential Equations, 3 (2016) 205-216.


\bibitem{darania-mara} P. Darania,  multi-step collocation method for
nonlinear delay integral equations, Sahand Communications in
Mathematical Analysis (SCMA), 3 (2016) 47-65.

\bibitem{darania-pishbin} P. Darania and S. Pishbin, Stability Analysis of the High-order multi-step Collocation Method for the Functional
Integral Equations with Constant Delays, Bol. Soc. Paran. Mat. (3s.) v. 2022 (40), 1–11.


\bibitem{lin} M.Q. Xu, Y.Z. Lin and Y.H.Wang, A new algorithm for nonlinear fourth
order multi-point boundary value problems, Applied Mathematics and
Computation, 274: 163-168, 2016.

% ----------------------------------------------------------------------


 \bibitem{ali-brunner} I. Ali, H. Brunner and T. Tang, Spectral
methods for pantograph-type differential and integral equations with
multiple delays, Front. Math. China, 4 (2009) 49-61.
%----------------------------------------------------------------------
\bibitem{baker} C.T.H. Baker, M.S. Derakhshan, Convergence and
stability of quadrature methods applied to Volterra equations with
delay, IMA J. Numer. Anal. 13 (1993) 67-91.
%% ----------------------------------------------------------------------
\bibitem{bookbrunner} H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations,
 Cambridge University Press, 2004.
%---------------------------------------------------------------------------------
\bibitem{brunner} H. Brunner, High-order collocation methods for singular Volterra functional equations of neutral
 typr, Applied Numerical Mathematics, 57 (2007) 533-548.
% ----------------------------------------------------------------------
\bibitem{brunneriter} H. Brunner, Iterated collocation methods for  Volterra integral equations with
delay arguments, Math. Comp. 62 (1994) 581-599.
%----------------------------------------------------------------------
\bibitem{brunner1} H. Brunner, The numerical solution of weakly singular Volterra functional
 integro-differential equations with variable delays, Comm. Pure Appl. Anal. 5 (2006) 261-276.
% ----------------------------------------------------------------------
\bibitem{khasi} M. Khasi, F. Ghoreishi and M. Hadizadeh, Numerical analysis of a high order
method for state-dependent delay integral equations, Numer.
Algorithms, 66  (2014) 177-201.
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

 \end{thebibliography}

\end{document}

