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

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2014 (2014), No. 54, pp. 1--9.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2014 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2014/54\hfil Block-pulse functions]
{Block-pulse functions and their applications to solving systems of
 higher-order nonlinear Volterra integro-differential equations}

\author[A. Ebadian, A. A. Khajehnasiri \hfil EJDE-2014/54\hfilneg]
{Ali Ebadian, Amir Ahmad Khajehnasiri}  % in alphabetical order

\address{Ali Ebadian \newline
Department of Mathematis, Urmia University, Urmia, Iran}
\email{a.ebadian@urmia.ac.ir}

\address{Amir Ahmad Khajehnasiri \newline
Department of Mathematis, Urmia University, Urmia, Iran. \newline
Department of Mathematics, Payame Noor University, 
PO Box 19395-3697 Tehran, Iran}
\email{a.khajehnasiri@gmail.com}

\thanks{Submitted  May 4, 2013. Published February 21, 2014.}
\subjclass[2000]{45G10, 45D05}
\keywords{Operational matrix; 
  Volterra integral equations; \hfill\break\indent block-pulse functions}

\begin{abstract}
 The operational block-pulse functions, a well-known method for solving functional
 equations, is employed to solve a system of nonlinear Volterra integro-differential
 equations. First, we present the block-pulse operational matrix of integration,
 then by using these matrices, the nonlinear  Volterra high-order
 integro-differential equation is reduced to an algebraic system.
 The benefits of this method is low cost of setting up the equations without
 applying any projection method such as Galerkin, collocation, etc.
 The results reveal that the method is very effective and convenient.
\end{abstract}

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


\section{Introduction}

Systems of integro-differential equations are a well-known mathematical
tool for representing physical problems. Historically, they have achieved
great popularity among the mathematicians and physicists in formulating
boundary value problems of gravitation, electrostatics, fluid dynamics
and scatering.

The concept of the block-pulse functions (BPFs) was first introduced to
electrical engineers by Harmuth. Then several researchers
(Gopalsami and Deekshatulu, 1997 \cite{ref8}, Sannuti, 1977 \cite{ref9},
 Riad,  1992 \cite{ref402}, Chen and Tsay, 1977 \cite{ref10})
discussed the BPFs and their operational matrix \cite{ref11, ref12}.

The aim of this work is to present a numerical method for approximating
the following system of   nonlinear Volterra  integro-differential equations
of  order $r$ $(r \geq 1)$:
\begin{equation}\label{f1}
u_i^r(x)+\sum_{k=1}^r B_{r-k}u_i^{(r-k)}(x)
=g_i(x)+\lambda\int_0^xk_i(x,t,F(u_i(x)))dt, \quad i=1,2,\dots,n,
\end{equation}
with initial conditions
\begin{equation}
u_i(0)=u_i,
\end{equation}
where the parameters $\lambda$ and functions $g_i(x)$, $k_i(x,t,F(u_i(t))$
are known and belong to $L^{2}[0,1)$, and  $u(x)$ is an unknown function
and $B_i$ $(i=0,1,\dots,n)$  are $n\times n$ matrices.
In this work, we consider that, the nonlinear function has the  form
\[
F(u_i(t))=(u_i(t))^p ,
\]
where $p$ is a positive integer.

Systems of integro-differential equations have a major role in the fields
of science and engineering \cite{ref1201, ref00, ref000, ref0001}.
The initial value problem for a nonlinear system of integro-differential
equations was used to model the competition between tumor cells and the
immune system \cite{ref121}. There are various techniques for solving a
system of integral or integro-differential equation, e.g. Galerkin
 method \cite{ref}, Adomian decomposition method (ADM)
 \cite{ref1, ref2}, rationalized Haar functions method \cite{ref3} and
variational iteration method (VIM) \cite{ref4}, He's homotopy perturbation
method (HPM) \cite{ref5,ref6}, Tau method \cite{ref7}, differential
transform method \cite{ref1111}, and Maleknejad in \cite{ref0000}
 have applied  Bernstein operational matrix  for solving a system of
high order linear Volterra-Fredholm integro-differential equations, etc.

This article is organized as follows:
In Section 2, we introduce  BPFs and their properties.
In Section 3, the operational matrix of integration is derived.
Section 4 introduces  Application of the method. Some numerical results has
been presented in section 5 to show accuracy and advantage of the
proposed method. Finally, some concluding remarks are given in section 6.

\section{Properties of  block-pulse functions}

An $m$-set of  BPFs is defined as follows:
\begin{equation}\label{00}
\phi_i(t) = \begin{cases}
1, & ih \leqslant t < (i+1)h, \\
0, & \text{otherwise},
\end{cases}
\end{equation}
where $i=1,2,\dots,m-1$ with positive integer values for $m$, and
$h=T/m$, and $m$ are arbitrary positive integers.
There are some properties for BPFs, e.g. disjointness, orthogonality,
 and completeness.
\smallskip

\noindent\textbf{Disjointness.}
The  block-pulse functions are disjoint with each other; i.e.,
\begin{equation}\label{0000}
\phi_i(t)\phi_{j}(t)= \begin{cases}
\phi_i(t), & i=j, \\
0 ,& i\neq j,
\end{cases}
\end{equation}
where $i,j=0,\dots,m-1$.
\smallskip

\noindent\textbf{Orthogonality.}
The  block-pulse functions are orthogonal with each other; i.e.,
\begin{equation}
\int_0^{T}\phi_i(t)\phi_{j}(t) dt =\begin{cases}
h, & i=j,\\
0 ,& \text{otherwise},
\end{cases}
\end{equation}
in the region of $ t\in[0,T)$,  where $i,j=1,2,\dots, m-1$. This property
is obtained from the disjointness property.
\smallskip

\noindent\textbf{Completeness.}
For every $f\in L^{2}([0,1))$  when $m$  go to infinity, Parseval identity holds:
\begin{equation}
\int_0^{1}f^{2}(t) dt = \sum_{i=0}^\infty  f_i^{2} \|\phi_i(t) \|^{2},
\end{equation}
 where
\[
 f_i=\frac{1}{h}\int_0^{1} f(t)\phi_i (t)dt.
\]
\medskip

The set of BPFs may be written as a $m$-vector $\Phi(t):$
\begin{equation}\label{6666}
\Phi(t)=[\phi_0(t),\dots,\phi_{m-1}(t)]^{T},
\end{equation}
where $ t\in [0,1)$.
From the above representation and disjointness property, it follows:
\begin{gather*}
 \Phi(t)\Phi^{T}(t) = \begin{bmatrix}
\phi_0(t) & 0 & \dots & 0 \\
0 & \phi_{1}(t) & \dots & 0  \\
\vdots & \vdots & \ddots & \vdots\\
0&0& \dots & \phi_{m-1}(t)\\
\end{bmatrix} \\
 \Phi^{T}(t)\Phi(t)=1,\\
\Phi(t)\Phi(t)^{T}V=\tilde{V}\Phi(t),
\end{gather*}
where $V$ is an $ m$-vector and $V=\operatorname{diag}(V)$.
Moreover, it can be clearly concluded that for every $m\times m$ matrix $A$:
\begin{equation}\label{1011}
\Phi^{T}(t)A\Phi(t)=\widehat{A^{T}}\Phi(t),
\end{equation}
where $A$ is an $ m$-vector with elements equal to the diagonal entries of matrix
$A$.

\subsection{Functions approximation}

A function $f(t)\in L^2([0,1))$ may be expanded by the BPFs as:
\begin{equation}\label{2.1}
f(t)\simeq \sum_{i=0}^{m-1}  f_{i_1}\phi_i(t)=F^{T}\Phi(t)=\Phi^{T}(t)F,
\end{equation}
where $F$ is a $m$-vector given by
\begin{gather}
F=[f_0,\dots,f_{m-1}]^{T}, \\
\label{2.6}
\Phi (t)=[\phi_{1}(t),\phi_{2}(t),\dots ,\phi_{m-1}(t)]^{T},
\end{gather}
the block-pulse coefficients $f_i$ are obtained as
\begin{equation}\label{2.2}
f_i=\frac{1}{h}\int_{ih}^{(i+1) h} f(t)dt,
\end{equation}
such that error between $f(t)$, and its block-pulse expansion \eqref{2.1}
in the region of $t \in [0,1)$
\begin{equation}
\varepsilon=\int_0^{1}\Big(f-\sum_{i=0}^{m-1}f_i \phi_i(t)\Big)^2dt,
\end{equation}
is minimal.
Now assume  $K(x,t) \in  L^{2}([0,1)\times [0,1))$ may
be approximated with respect to BPFs such as:
\begin{equation}
k(x,t)=\Phi^{T}(x)K \Phi(t),
\end{equation}
where $\Phi(x)$ and $\Phi(t)$ are BPFs vectors of dimension $m_1$ and $m_2$,
respectively, and K is a $m_1\times m_2 $ one dimensional block-pulse
coefficients matrix with $k_{ij}$, $i=0,\dots ,m_{1}-1$, $j=0,\dots,m_{2}-1$
as  follows:
\begin{equation}
k_{ij}=m_{1}m_{2}\int_0^{1}\int_0^{1}k(x,t)\phi_i(x)\phi_{j}(t)\,dx\,dt.
\end{equation}
Also, the positive integer powers of a function $f(s)$ may be approximated by
BPFs as
\[
[u(t)]^{p}= [\Phi^T(t) U]^{p}=\Phi^T(t) \Lambda,
\]
where $\Lambda$  is a column vector, whose elements are $p$th power of the
elements of the vector $U$.

\subsection{Block-pulse functions series}

The function $x^k$, $x \in [0,1)$, $k \in \mathbb{N}$ can be
approximated as a BPF series of size $m$. Indeed, from \eqref{2.1}
and \eqref{2.2}, we have
\begin{equation}
x^k\simeq \sum_{i=0}^{m-1} f_{k}(i)\phi(x),
\end{equation}
where
\begin{equation}
f_{k}(i)=\frac{1}{h}\int_{ih}^{(i+1) h} t^k dt=\frac{1}{h(k+1)}
[ ((i+1)h)^{k+1}-(ih)^{k+1} ].
\end{equation}
Therefore,
\begin{equation}
x^k\simeq \frac{1}{h(k+1)}\sum_{i=0}^{m-1}[ ((i+1)h)^{k+1}-(ih)^{k+1}]\Phi _i(x),
\end{equation}
and in matrix form
\begin{equation}\label{2.4}
x^k\simeq \frac{1}{h(k+1)}Y^{T}_{k}\Phi _{m}(x),
\end{equation}
where
\[
Y_{k}^{T}=\sum_{i=0}^{m-1}[ ((i+1)h)^{k+1}-(ih)^{k+1}].
\]

\section{Operational matrix of integration}

We compute $\int_0^{t} \Phi_id\tau  $ as
\begin{equation}\label{1111}
\int_0^{t}\Phi _i(\tau)d \tau =\begin{cases}
0,& t \leq ih,\\
t-ih & ih\leq t< (i+1)h, \\
h & (i+1)h\leq t<1.
\end{cases}
\end{equation}
Then \eqref{1111} can be written as
\begin{equation}\label{4444}
\int_0^{t}\Phi _i(\tau)d \tau=(t-ih)\Phi _i(t)+h\sum_{j=i+1}^{m-1}\Phi _{j}(t).
\end{equation}
From \eqref{2.4} we have
\begin{equation}\label{2222}
x \simeq \frac{1}{2h} \sum_{i=0}^{m-1}[   ((i+1)h)^{2}-(ih)^{2} ] \Phi _i(t).
\end{equation}
Substituting \eqref{2222} and \eqref{0000} into \eqref{4444}, and by using
 orthogonal property, for $ 0\leq i < m $, we have
\begin{align*}
 \int_0^{t}\Phi _i(\tau)d \tau
 &= \frac{1}{2h}\sum _{j=0}^{n-1}\left[ ((j+1)h)^{2}-(jh)^{2}\right]
 \Phi_{j}(t)\Phi_i(t)-ih\Phi_i(t)+h\sum_{j=i+1}^{m-1}\Phi _{j}(t)\\
 &= \frac{1}{2h}[((i+1)h)^{2}-(ih)^{2}]\Phi_i(t)-ih\Phi_i(t)
 +h\sum_{j=i+1}^{m-1}\Phi_{j}(t)\\
 &= \frac{h}{2}\Phi _i+h\sum_{j=i+1}^{m-1}\Phi_{j}(t).
\end{align*}
The integration of the vector $\Phi(t)$ defined in \eqref{6666} may be
obtained as
\begin{equation}\label{2.3}
\int_0^{t}\Phi(\tau)d \tau  \simeq \Upsilon \Phi(t),
\end{equation}
where $\Upsilon$ is called operational matrix of integration
 which can be represented by
\[
\Upsilon=\frac{h}{2}
\begin{pmatrix}
1 & 2 & 2 &\dots & 2\\
0 & 1 & 2 &\dots & 2 \\
\vdots & \vdots & \vdots &\ddots &\vdots\\
0 & 0 & 0 & \dots &1
\end{pmatrix},
\]
and their integrals in the matrix form
\[
\begin{pmatrix}
\int \Phi_0 \\
\int \Phi_{1} \\
\vdots \\
\int \Phi_{m-1}
\end{pmatrix}
\simeq \frac{h}{2} \begin{pmatrix}
1 & 2 & 2 &\dots & 2\\
0 & 1 & 2 &\dots & 2 \\
\vdots & \vdots & \vdots &\ddots &\vdots\\
0 & 0 & 0 & \dots &1
\end{pmatrix}
\begin{pmatrix}
 \Phi_0 \\
 \Phi_{1} \\
\vdots \\
 \Phi_{m-1}
\end{pmatrix},
\]
or in more compact form
\begin{equation}\label{7777}
\int_0^{t}\Phi_{m}(\tau)d \tau  \simeq \Upsilon \Phi_{m}(t),
\end{equation}
By using this matrix, we can express the integral of a function $f(t)$
into its block pulse series
\begin{equation}
\int_0^{t} f_{m}(\tau)d \tau \simeq \int_0^{t}  F^{T}\Phi_{m}
(\tau)d \tau\simeq F^{T} \Upsilon\Phi_{m}(t).
\end{equation}

\section{Application of the method}

In this section, we calculate $U_i^{k-r}(x)$ by using
\begin{equation}\label{3.1}
 u_i^r(x)=\sum_{i=0}^{m-1}u_i\Phi(x)=U_i^{T}\Phi_{m}(x).
 \end{equation}
 Now integrating   from $0$ to $t$ and using \eqref{7777} we obtain
 \begin{equation}
U_i^{r-1}(x)=U\Upsilon \Phi_{m}(x)+U_0^{r-1}(x)
\end{equation}
The $k-th $ integration of \eqref{3.1} yields
\begin{equation}\label{3.3}
U_i^{r-k}(x)=U_i\Upsilon ^k\phi_{m}(x)
+\sum_{i=1}^kU_0^{r-i}\frac{t^{k-i}}{(k-i)!} \quad k=1,2,\dots ,r.
\end{equation}
From \eqref{2.4} we have
\begin{equation}\label{3.4}
\frac{t^{k-i}}{(k-i)!}\simeq \frac{1}{h(k-i+1)!}Y_{k-i}^{T}\Phi_{m}(x)
\end{equation}
Substituting \eqref{3.4} in \eqref{3.3} we obtain
\begin{equation}\label{1010}
U_i^{r-k}(x)=U_i\Upsilon ^k\Phi_{m}(t)+Z_{k}\Phi_{m}(x),
\end{equation}
where
\begin{equation}
Z_{k}=\frac{1}{h}\sum_{i=1}^k\frac{1}{h(k-i+1)!}U_0^{r-i}Y_{k-i}^{T}
\end{equation}
is an $n \times m $ constant matrix.

Now, we solve  the system of nonlinear  Volterra high-order integro-differential 
equations by using BPFs. As we show before, we can write
\begin{equation}\label{8888}
\begin{gathered}
g_i(x)=G_i^{T}\Phi_{m}(x),\\
 u_i^r(x)=U_i^{T}\Phi_{m}(x),\\
[u_i(x)]^{p} =  \Phi_{m}^T(x) \Lambda,\\
k(x,t)=\Phi^{T}(x)K \Phi(t),
\end{gathered}
\end{equation}
where the $m$-vectors $U, G, \Lambda, $ and matrix $K$ are BPFs 
coefficients of $u(x), g(x),[u(t)]^{p}, $ and $K(x,t)$ respectively, 
$\Lambda $ is a column vector whose elements are $p$th power of the 
elements of the vector U.
To approximate the integral equation \eqref{f1}, from \eqref{8888} and \eqref{1010}
 we get
\[
u_i^r(x)+\sum_{k=1}^r B_{r-k}u_i^{(r-k)}(x)=g_i(x)+\int_0^xk_i(x,t,F(u_i(x)))dt \quad
i=1,2,\dots,n.
\]
Now the second part of equation
\begin{align*}
U_i^{T}\Phi _{m} (x)+\sum_{k=1}^rB_{r-k}(U_i\Upsilon ^k+Z_{k})\Phi _{m}(x)
 &= G_i^{T}\Phi_{m} (x)+\Phi_{m}^{T} (x)K\int_0^x\Phi(t)\Phi^{T}(t)\Lambda\\
 &=   G_i^{T}\Phi _{m}(x)+\Phi_{m}^{T} (x)K\tilde{\Lambda}\int_0^x\Phi_{m}(t)dt\\
 &=  G_i^{T}\Phi_{m} (x)+\Phi_{m}^{T} (x)K\tilde{\Lambda}\Upsilon\Phi_{m}(x).
\end{align*}
If we put $ A=K\tilde{\Lambda}\Upsilon$ then it can be written from \eqref{1011},
\begin{align*}
U_i^{T}\Phi_{m}(x)+\sum_{k=1}^rB_{r-k}u_i\Upsilon ^k\Phi_{m}(x)=G_i^{T}\Phi_{m} (x)+\widehat{A^{T}}\Phi_{m}(x),
\end{align*}
hence, we have
\begin{equation}
U_i^{T}+\sum_{k=1}^rB_{r-k}u_i\Upsilon ^k=G_i^{T}+\widehat{A^{T}}.
\end{equation}
It can be written as:
\begin{equation}\label{c92}
AU=F
\end{equation}
where $A$ and $F$ are the combination of block-pulse coefficient matrix and 
$U$ can be obtained from Newton-Raphson method for solving nonlinear systems.

\section{Numerical examples}

To illustrate the effectiveness of the proposed method in the present paper, 
several test examples are carried out in this section.

\begin{example} \label{examp1} \rm
Consider the  nonlinear  Volterra  integro-differential equations  
problem with initial conditions \cite{ref5},
\begin{equation}
\begin{gathered}
u'(x)-1+\frac{1}{2}{{v}^{' }}^{2}(x)= \int_0^x((x-t)v(t)+v(t)u(t))dt, \\
v'(x)-2x=\int_0^x((x-t)u(t)-v^{2}(t)+u^{2}(t))dt, \\
u (0)=0,\quad v(0)=1.
\end{gathered}
\end{equation}
\end{example}

The exact solutions are $u(x) =\sinh(x)$,  $v(x) =\cosh(x)$. 
The numerical results  obtained with BPFs  are presented  in Table \ref{Tab1}  
and Figure \ref{shir}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig1}
\end{center}
\caption{Comparison of the exact solution  and the present method}
\label{shir}
\end{figure}

\begin{example} \label{examp2}\rm
consider the system of two nonlinear integro-differential 
equations with boundary conditions \cite{ref5},
\begin{equation}
\begin{gathered}
u''(x)=1-\frac{1}{3}x^{3}-\frac{1}{2}{v'}^{2} (x)
 +\frac{1}{2}\int_0^x (u^{2}(t)+v^{2}(t))dt,\\
v''(x)=-1+x^{2}-xu(x)+\frac{1}{4}\int_0^x(u^{2}(t)-v^{2}(t))dt, \\
u(0)=1,\quad u'(0)=2,\quad v(0)=-1,\quad  v'(0)=0.
\end{gathered}
\end{equation}
\end{example}

\begin{table}[htb]
\caption{Numerical results of Example \ref{examp1}}\label{Tab1}
\begin{center} \vskip -10pt
{\scriptsize
\begin{tabular}{ccccc}\hline
$x$ &$(u_{\rm exact}(x), v_{\rm exact}(x))$ & $m=8$ & $m=16$ & $m=32$     \\
\hline
$0.0$ & $(0.00000, 1.00000)$ & $(0.00145, 0.88425)$ & $(0.00115, 0.91196)$
 & $(0.00025, 1.00056)$    \\
$0.1$ & $(0.10016, 1.00500)$ & $(0.10727, 1.01570)$ & $(0.101852, 1.009)$
 & $(0.10016, 1.00501)$ \\
$0.3$ & $(0.30452, 1.04533)$ & $(0.38440, 1.00521)$
 & $(0.58562, 1.03122)$ & $(0.30450, 1.04533) $   \\
$0.5$ & $(0.52109, 1.12762)$ & $(0.50695, 1.108423)$ & $(0.52012, 1.12521)$ & (0.52108, 1.12762)
  \\
$0.7$ & $(0.75858, 1.25516)$ & $(0.74932, 1.20390)$ & $(0.75125, 1.25501)$
 & $(0.75857, 1.25510)$         \\
$0.9$ & $(1.02651, 1.43308)$ & $(1.09032, 0.30390)$ & $(1.02541, 1.43321)$
 & $(1.02651, 1.43308)$  \\ \hline
\end{tabular}
}\end{center}
\end{table}

The exact solutions are  $ u(x) =x+ e^x$ and  $v(x) =x- e^x$. 
Numerical results for this solution is presented  in Table \ref{Tab2}.

\begin{table}[ht]
\caption{Numerical results of Example \ref{examp2}}\label{Tab2}
\begin{center} \vskip -10pt
{\scriptsize
\begin{tabular}{ccccc}\hline
$x$ & $(u_{exact}(x), v_{exact}(x))$ & $m=8$ & $m=16$ & $m=32$   \\ \hline
$0.0$ & $(1, -1)$ & $(0.93915, -0.88365)$ & $(0.96251, -0.89936)$ & 
$(0.99251, -0.97936)$ \\
$0.1$ & $(1.02051, -1.00517)$ & $(1.01107, -1.02517)$ & $(1.02012, -1.01701)$
  &$(1.02052, -1.00701)$     \\
$0.3$ & $(1.64985, -1.04985)$ & $(1.61852, -1.00980)$ & $(1.64212, -1.04914)$
 & $(1.64812, -1.04984)$ \\
$0.5$ & $(2.14872,-1.14872)$ & $(2.12211,-1.14582)$ & $(2.14121, -1.14705)$ 
 & $(2.14821, -1.14725)$    \\
$0.7$ & $(2.71375,-1.31375)$ & $(2.71175,-1.31075)$ & $(2.71301,-1.31221)$ & 
 $(2.71371,-1.31321)$  \\
$0.9$ & $(3.35960,-1.55960)$ & $(3.35255,-1.54696)$ & $(3.35666,-1.55009)$ & 
$(3.35966,-1.55909)$    \\
      \hline
\end{tabular}
}\end{center}
\end{table}

\subsection*{Conclusion}
In this article, we  approximated the solution  of nonlinear Volterra 
integro-differential equations. To this end, we  used some orthogonal 
functions called Block-Pulse Functions. Finally, numerical examples reveal 
that the present method is very accurate and convenient for solving systems 
of high order linear and  nonlinear Volterra  integro-differential equations. 
The benefits of this method is low cost of setting up the equations without 
applying any projection method such as Galerkin, collocation, etc.
 Also, the linear system \eqref{c92} has a regular form which can help us 
for solving it.

\begin{thebibliography}{10}

\bibitem{ref1111} A. Arikoglu, I. Ozkol;
 \emph{Solutions of integral and integro-differential equation systems by using 
differential transform method}, Computers and Mathematics with Applications, 
\textbf{56} (2008) 2411-2417.

\bibitem{ref7} S. Abbasbandy, A. Taati;
 \emph{Numerical solution of the system of nonlinear Volterra integro-differential 
equations with nonlinear differential part by the operational Tau method 
and error estimation}, Journal of Computational and Applied Mathematics, 
\textbf{231} (2009) 106-113.

\bibitem{ref1201} M. I. Berenguer, A. I. Garralda-Guillem, M. Ruiz Galan;
\emph{An approximation method for solving systems of Volterra
integro-differential equations,} Applied Numerical Mathematics, 
\textbf{67} (2013) 126-135.

\bibitem{ref121} N. Bellomo, B. Firmani, L. Guerri;
\emph{Bifurcation analysis for a nonlinear system of integro-differential 
equations modelling tumor–immune cells competition,} 
Appl. Math. Lett, \textbf{12} (1999) 39-44.

\bibitem{ref2} J. Biazar, E. Babolian, R. Islam;
\emph{Solution of a system of Volterra integral equations of the first 
kind by Adomian method}, Appl. Math. Comput, \textbf{139} (2003) 249-258.

\bibitem{ref5} J. Biazar, H. Ghazvini, M. Eslami;
 \emph{He's homotopy perturbation method for systems of integro-differential 
equations}, Chaos Solitons. Fractals, \textbf{39} (2009) 1253-1258.

\bibitem{ref10} C. F. Chen, Y. T. Tsay, T. T. Wu;
\emph{Walsh operational matrices for fractional calculus and their 
 application to distributed systems,} Journal of the Franklin Institute, 
\textbf{303} (1977) 267-284.

\bibitem{ref11} A. Ganti, Prasada Rao;
 \emph{Piecewise Constant Orthogonal Functions and their Application to System
 and Contro,} Springer-Verlag, 1983.

\bibitem{ref8} G. Gopalsami, B. L. Deekshatulu;
 \emph{Comments on Design of piecewise stant gains for optimal
control via Walsh functions,} IEEE Trans. on Automatic Control, 
\textbf{123} (1997) 461-462.

\bibitem{ref00} M. A. Jaswon, G. T. Symm;
\emph{Integral Equation Methods in Potential Theory and Elastostatics,}
Academic Press, London, 1977.

\bibitem{ref000} A. Kyselka;
 \emph{Properties of systems of integro-differential equations in the statistics 
of polymer chains,} Polym. Sci. USSR, \textbf{19} (11) (1977) 2852-2858.

\bibitem{ref12} B. M. Mohan, K. B. Datta;
\emph{Orthogonal Function in Systems and Control,} 1995.

\bibitem{ref} K. Maleknejad, M. Tavassoli Kajani;
 \emph{Solving linear integro-differential equation system by Galerkin methods 
with hybrid functions}, Appl. Math. Comput, \textbf{159} (2004) 603-612.

\bibitem{ref3} K. Maleknejad, F. Mirzaee, S. Abbasbandy;
 \emph{Solving linear integro-differential equations system by using 
rationalized Haar functions method,} Appl. Math. Comput, 
\textbf{155} (2004) 317-322.

\bibitem{ref0000} K. Maleknejad,  B. Basirat, E. Hashemizadeh;
\emph{A Bernstein operational matrix approach for solving a system of high
order linear Volterra–Fredholm integro-differential equations}, 
Mathematical and Computer Modelling, \textbf{55} (2012) 1363-1372.

\bibitem{ref402} R. Riad;
 \emph{Solution of system of high-order differential equation with constant 
coefficients via Block-pulse functions}, Annales univ Sci, 
\textbf{13} (1992) 11-20.


\bibitem{ref9} P. Sannuti;
 \emph{Analysis and synthesis of dynamic systems via Block Pulse functions, 
IEEE Proceeding,} \textbf{124} (1977) 569-571.

\bibitem{ref1} H. Sadeghi Goghary, Sh. Javadi, E. Babolian;
 \emph{Restarted Adomian method for system of nonlinear Volterra integral equations},
 Appl. Math. Comput, \textbf{161} (2005) 745-751.

\bibitem{ref0001} P. Schiavane, C. Constanda, A. Mioduchowski;
 \emph{Integral Methods in Science and Engineering, Birk\"auser,} Boston, 2002.

\bibitem{ref4} S. Q. Wang, J. H. He;
\emph{Variational iteration method for solving integro-differential equations},
Phys. Lett. A, \textbf{367} (2007) 188-191.

\bibitem{ref6} E. Yusufoglu (Agadjanov);
 \emph{An efficient algorithm for solving integro-differential equations system},
 Appl. Math. Comput, \textbf{192} (2007) 51-55.

\end{thebibliography}
\end{document}