\documentclass[twoside]{article}
%\usepackage{graphics}
\usepackage{graphicx}
\usepackage{latexsym}
\usepackage{amssymb}
\usepackage{amsmath}
%\usepackage{cite}
\usepackage{BSPMsample}
\begin{document}
\title[Numerical  Simulations  of Two-Species Reaction-Diffusion System]
{Numerical  Simulations of Two Components  Ecological
Models with Sinc Collocation Method}
%between [] goes the title that will appear on the top of every
%odd page, between {} goes the real title.
%\thanks{...} put your grant

\author{A. Barati \thanks{Corresponding author}  and A. Atabaigi}
%between [] goes the authors that will appear on the top of every
%even page, between {} goes the real author.

\address{Ali Barati\\Eslam Abad Gharb Faculty of Engineering \\
, Razi University, Kermanshah \\ Iran\\
\email: alibarati@razi.ac.ir}
\address{Ali Atabaigi\\Department of Mathematics, Faculty of Science\\ Razi University, Kermanshah \\ Iran\\
\email: a.atabeigi@razi.ac.ir}
\maketitle

\begin{abstract}
It is well known that reaction-diffusion systems can describe the interactions of different species in ecological systems. In this work, we propose a numerical method  for solving a
general two-species reaction-diffusion system. The method comprises a standard
finite difference to discretize  in time direction and
Sinc collocation method in spatial direction. A
series of numerical experiments demonstrate  the accuracy and good performance of the algorithm. The biological significance of the numerical results was also discussed and plotted in figures.
\end{abstract}
%%put here Key words

\keywords Reaction-diffusion systems, Predator-prey, Sinc collocation method.
\tableofcontents
%% put here the subject class
\2000mathclass{65M22, 65N30, 35F20}


\section{Introduction}\label{sec1}
Reaction-diffusion equations may be used to model many natural
phenomena and practical problems which arise from a variety of
disciplines such as physics, chemistry, medicine and so on. These
equations are a useful tool for modeling the spatiotemporal
dynamics of populations in ecology \cite{can,hal}, The study of
reaction-diffusion problems in ecological context have gained a
huge amount of scientific interest, due to their practical
relevance and emergence of some interesting phenomena such as
spatial patterns, oscillating solutions, phase planes and multiple
steady states to mention a few.\\
In recent years, there have been
considerable interests in spatial and temporal behavior of
interacting species in ecosystems.
 Pattern formation study in
reaction-diffusion systems is a very active research area.
 Since
Turing \cite{Tur} first proposed reaction-diffusion theory to
describe the range of spatial patterns observed in the developing
embryo, reaction-diffusion systems have been studied extensively
to explain patterns in fish skin, mammalian coat markings,
phyllotaxis, predator-prey systems, terrestrial vegetation,
plankton, intertidal communities and so on (see
Ref. \cite{duf,sun,sun1,liu,Wan}). \\
Predator-prey model is one
of the most important population dynamical models. There are many
factors which affect population dynamics in predator-prey models.
One crucial component of predator-prey relationships is
predator-prey interaction (also called functional response), which
can be classified in to many different types, such as Holling I-IV
types, Hassell-Varley type, Beddington-DeAngelis type,
Crowley-Martin type, and etc.\\
Predator-prey systems have been studied by many researchers in
various forms. For instance, in bacteria ecology, computer
simulations of complex spatiotemporal patterns \cite{baek,garv} of
Bacillus subtilis based on deterministic models \cite{Mi}, Allee
effect of patchy invasion on predator-prey dynamics \cite{pet}.
 Jiang et al in \cite{jia} consider a predator-prey model
with Beddington-DeAngelis functional response subject to the
homogeneous Neumann boundary condition. Moreover, Wang et al.
\cite{wan} investigated the spatial pattern formation of a
predator-prey system with prey-dependent functional response of
Ivlev type and reaction-diffusion whereas the analysis of
predator-prey systems showing the Holling type II functional
response is examined in Garvie and Trenchea \cite{garv1}.
Moreover, we can point to other efficient numerical methods for
solving the reaction-diffusion equations arising in biology such
as \cite{kol,fat,mar,kol1,mar1}.\\
In this paper, we consider the following  general two-species
reaction-diffusion system
\begin{equation}\label{conEq}
\left\{\begin{array}{c} \frac{\partial u}{\partial t}= d_u
\frac{\partial ^{2}u}{\partial x ^{2}}+f(u,v), \hskip 1cm (x,t)\in
\Omega \times (0,\infty),
\\\\
\frac{\partial v}{\partial t}= d_v \frac{\partial ^{2}v}{\partial
x ^{2}}+g(u,v), \hskip 1cm (x,t)\in \Omega \times (0,\infty),
                                                 \end{array}
                                                 \right.
\end{equation}
subject to the initial conditions:
\begin{equation}\label{cond1}
u(x,0)=u_{0}(x), \hskip 1cm v(x,0)=v_0(x) \hskip 1cm  x \in
\overline{\Omega},
\end{equation}
and zero-flux boundary conditions
\begin{equation}\label{cond2}
\frac{\partial u}{\partial \nu }=\frac{\partial v}{\partial \nu
}=0,\hskip .75cm (x,t) \in
\partial \Omega \times (0,\infty)  ,
\end{equation}
where u(x,t) and v(x,t) represent the species densities, $\Omega$
is a bounded  region with the smooth boundary $\partial \Omega$ in
$\mathbb{R} $ ( here we assume $\Omega=(a,b)$), $\nu$ is the outward unit normal vector on
$\partial \Omega$. The positive constants $d_u$ and $d_v$ are the
 diffusion coefficients corresponding to $u$ and $v$, also $f$ and $g$ are reaction kinetics.\\
The homogeneous Neumann boundary condition means that model
(\ref{conEq}) is self-contained and has no population flux across
the boundary $\partial \Omega $.\\
We assume that the point $(u_s,v_s)$ is a positive equilibrium point  of
the homogeneous system
\begin{equation}
\frac{d u}{d t}=f(u,v),\hskip 1cm
 \frac{d u}{d t}=g(u,v),
\end{equation}
that is $f(u_s,v_s)=0, \hskip .25cm g(u_s,v_s)=0$\\
In \cite{Tur}, Turing concluded that the reaction-diffusion model (\ref{conEq}) may exhibit spatial patterns under the following two conditions: the equilibrium $(u_s,v_s)$ is linearly stable in the absence of diffusion; and the equilibrium becomes linearly unstable in the presence of diffusion. Such an instability is called a \textit{Turing instability} or \textit{diffusion-driven instability}.\\
  \indent In this paper, a numerical method applied  for solving  couple system (\ref{conEq}) which is based on Sinc basis
functions.  Sinc methods have been studied extensively and found
to be a very effective technique, particularly  for problems
  with singular solutions and those on unbounded domain.
  In addition, Sinc function seems that be useful for problems  that
  their solutions have  oscillatory behavior in domain.
  Sinc method originally introduced by Stenger \cite{st}
 which is based on the Whittaker-Shannon-Kotel' nikov sampling theorem for entire
functions. The books \cite{lu}  and \cite{st1} provide excellent
overviews of the existing Sinc methods for solving ODEs and
PDEs.\\
\indent In recent years, a lot of attentions have been devoted to
the study of the Sinc method to investigate various scientific
models. The efficiency of the Sinc method has been formally proved
by many researchers  Bialecki \cite{bi}, Rashidinia et al.
\cite{RMA,rash},  EL-Gamel \cite{ga}, Okayama et al. \cite{to} and Saadatmandi and Dehghan
\cite{sa}. \\
The paper is organized as follows. In section 2, we review some
basic facts about the Sinc approximation. In section 3, we
discretized the reaction-diffusion system  in temporal variable by means of implicit Euler
method and then we applied the Sinc-collocation method  for
solving of the arising system of nonlinear ordinary differential
equations in each time level. Some numerical
examples will be presented in section 4, and at the end we
conclude implementation, application and efficiency of the proposed
scheme.

\section{Notation and background}\label{sec2}
The goal of this section is to recall notations and definitions of
the Sinc function and state some known theorems that are important
for this paper.\\
 The Sinc function is defined on
$-\infty<x<\infty$ by
$$Sinc(x)=\left\{\begin{array}{c}
                                                         \frac{sin(\pi x)}{\pi x}  , \hskip 0.5cm x\neq0 , \\
                                                           1, \hskip 1.25cm x=0 .
                                                  \end{array}
     \right.$$
For $h>0$ we will denote the Sinc basis functions by
$$ S(k, h)(x)=sinc(\frac{x-kh}{h}),\hskip 1cm k=0,\pm1,\pm2,\ldots $$
let $f$ be a function defined on $\mathbb{R}$ then for $h>0$ the
series
$$C(f,h)(x)=\sum_{k=-\infty}^{\infty} f(kh)S(k,h)(x),$$
is called the Whittaker cardinal expansion of  $f$ whenever this
series converges. The properties of Whittaker cardinal expansions
have been studied and are thoroughly surveyed in Stenger
\cite{st1}. These properties are derived in the infinite strip
$D_{d}$ of the complex plane where $d>0$
$$
D_{d}=\{\zeta=\xi+i\eta: |\eta|<d\leq\frac{\pi}{2}\}.$$
Approximations can be constructed for infinite, semi-finite, and
finite intervals. To construct approximation
on the interval $(a,b)$,we consider the conformal map\\
\begin{equation}\label{22}
\phi(z)=\ln(\frac{z-a}{z-b}),
\end{equation}\\
which maps the eye-shaped region
$$D_{E}=\{z=x+iy;|\arg(\frac{z-a}{z-b})|<d\leq\frac{\pi}{2}\},$$
onto the infinite strip $D_{d}$.\\
For the Sinc method, the basis functions on the interval $(a,b)$
for $ z\in D_{E}$  are derived from the composite translated Sinc
function:
\begin{equation}\label{22}
S_{k}(z)=S(k, h)\circ \phi(z)=sinc(\frac{\phi(z)-kh}{h}).
\end{equation}
The function
$$z=\phi^{-1}(\omega)=\frac{a+be^{\omega}}{1+e^{ \omega}},$$
is an inverse mapping of $\omega=\phi(z)$. We define the range of
$\phi^{-1}$ on the real line as
$$\Gamma=\{\psi (u)=\phi^{-1}(u) \in D_{E}: -\infty <u<\infty\}=(a,b).$$
The sinc grid points $z_{j}\in (a,b)$ in $D_{E}$ will be denoted
by $x_{j}$ because they are real. For the evenly spaced nodes
$\{jh\}_{j=-\infty}^{\infty}$ on the real line, the image which
corresponds to these nodes is denoted by
\begin{equation}\label{24}
 x_{j}=\phi^{-1}(jh)=\frac{a+be^{jh}}{1+e^{ jh}}, \hskip 1cm
 j=0,\pm1,\pm2,... .
\end{equation}
 \textbf{Definition 1.} Let $B(D_{E})$ is the class of functions
$f$
which are analytic in $D_{E}$ such that\\
\begin{equation}\label{def}
\int_{\psi(u+\Sigma)}|f(z)|dz\rightarrow 0, \hskip .75cm as \hskip
.25cm u \rightarrow\pm\infty
\end{equation}
where $\Sigma=\{i\eta: |\eta|<d\leq\frac{\pi}{2}\}$ and satisfy \\
\begin{equation}\label{def1}
N(f)\equiv \int_{\partial D_{E}}|f(z)|dz<\infty,
\end{equation}
where $\partial D_{E}$ represents the boundary of $D_{E}$.\\
\textbf{Definition 2.} Let $L_{\alpha}(D_{E})$ be the set of all
analytic function $u$ in $D_{E}$, for which there exists a
constant $C$ such that
\begin{equation}\label{def1}
|u(z)|\leq C \frac{|\rho (z)|^{\alpha}}{[1+|\rho (z)|
]^{2\alpha}}, \hskip .5 cm z\in D_{E},\hskip .25cm 0< \alpha \leq
1.
\end{equation}
where $\rho (z)=e^{\phi (z)}.$\\
 \textbf{Theorem 1.}(Stenger \cite{st1}) If $\phi^{'}u
\in B(D_{E})$,  and let
$$\sup_{\frac{-\pi}{h}\leq t\leq
\frac{\pi}{h}}\left|\left(\frac{d}{dx}\right)^{l}e^{it\phi(x)}\right|\leq
C_{1}h^{-l}, \hskip 1cm x \in \Gamma,$$ for $l=0,1,\ldots,m$ with
$C_{1}$ a constant depending only on $m$ and $\phi$. If $u \in
L_{\alpha}(D_{E})$ then taking $h=\sqrt{\pi d/\alpha N}$ it
follows that


$$\sup _{x \in
\Gamma}\left|u^{(l)}(x)-\left(\frac{d}{dx}\right)^{l}
\sum_{k=-N}^{N}u(x_{k})S_{k}(x)\right|\leq C N^{(l+1)/2}\exp(-(\pi
d\alpha N)^{1/2}),$$ where $C$ is a constant depending only on $u,
d, m, \phi$ and
$\alpha$.\\
 The Sinc-collocation method
requires that the derivatives of composite Sinc function be
evaluated at the nodes. We need to recall the following lemma.\\
\textbf{Lemma 1.}(Lund and Bowers \cite{lu}) Let $\phi$ be the
conformal one-to-one mapping of the simply connected domain
$D_{E}$ onto $D_{d}$, given by (\ref{22}). Then

\begin{equation}\label{25}
\delta_{kj}^{(0)}=[S(k,h)\circ
\phi(x)]|_{x=x_{j}}=\left\{\begin{array}{c}
                                                            1, \hskip 0.5cm k=j , \\
                                                           0, \hskip 0.5cm k\neq j ,
                                                  \end{array}
     \right.
 \end{equation}

 \begin{equation}\label{26}
\delta_{kj}^{(1)}=h\frac{d}{d\phi}[S(k,h)\circ
\phi(x)]|_{x=x_{j}}=\left \{ \begin{array}{c}
                                                                             0, \hskip 0.5cm k=j, \\
                                                                 \frac{(-1)^{(j-k)}}{j-k}, \hskip 0.5cm k\neq j ,
                                                                            \end{array}
\right.
\end{equation}

\begin{equation}\label{27}
\delta_{kj}^{(2)}=h^{2}\frac{d^{2}}{d\phi^{2}}[S(k,h)\circ
\phi(x)]|_{x=x_{j}}=\left \{ \begin{array}{c}
                                                                             \frac{-\pi^{2}}{3}, \hskip 0.5cm k=j, \\
                                                                 \frac{-2(-1)^{(j-k)}}{(j-k)^{2}}, \hskip 0.5cm k\neq j,
                                                                          \end{array}
\right.
\end{equation}
in relations (\ref{25}-\ref{27}) $h$ is step size and $x_{j}$ is
sinc grid given by (\ref{24}).
It is convenient to define the
following matrices:
\begin{equation}\label{28}
\textbf{I}^{(l)}=[\delta_{kj}^{(l)}], l=0,1,2,
\end{equation}
where $\delta_{kj}^{(l)}$ denotes the $(k,j)$th element of the
matrix $\textbf{I}^{(l)}$. Note that the matrix $\textbf{I}^{(2)}$
and $\textbf{I}^{(1)}$ are symmetric and
 skew-symmetric matrices respectively, also $\textbf{I}^{(0)}$ is identity
 matrix.



\section{Description of method}\label{sec3}
First, we discretize reaction-diffusion system (\ref{conEq}) in
time direction by means of the implicit Euler method with uniform
step size $\Delta t$ as \\
$$u_{t}=\frac{u^{n+1}-u^{n}}{\Delta t}, \hskip .75cm
v_{t}=\frac{v^{n+1}-v^{n}}{\Delta t},$$
 where $u^{n}=u(x,t_{n}), u^{0}=u(x,0), \hskip .25cm v^{n}=v(x,t_{n}), v^{0}=v(x,0)$ and $t_{n}=n\Delta t ,n=1,2,...$,

to get the following system of nonlinear ordinary differential
equations:
\begin{equation}\label{ode}
\left\{\begin{array}{c}  \frac{u^{n+1}-u^{n}}{\Delta t}= d_u
D^{2}u^{n+1}+f(u^{n+1},v^{n+1}) ,
\\\\
\frac{v^{n+1}-v^{n}}{\Delta t}= d_v
D^{2}v^{n+1}+g(u^{n+1},v^{n+1}),
                                                 \end{array}
                                                 \right.
\end{equation}
where $u^{n+1}= u(x,t_{n+1})$ and $v^{n+1}= v(x,t_{n+1})$ are the
solutions of Eqs.(\ref{ode}) at
  $(n+1)$th  time level and $D=\frac{d}{dx}$.\\
Now we can rewrite equation (\ref{ode}) in the following form
\begin{equation}\label{ode1}
\left\{\begin{array}{c} d_u
D^{2}\hat{u}+f(\hat{u},\hat{v})-\frac{1}{\Delta
t}\hat{u}=-\frac{u^{n}}{\Delta t} ,
\\\\
d_v D^{2}\hat{v}+g(\hat{u},\hat{v})-\frac{1}{\Delta
t}\hat{v}=-\frac{v^{n}}{\Delta t},
                                                 \end{array}
                                                 \right.
\end{equation}
where $\hat{u}=u^{n+1}, \hskip .25cm \hat{v}=v^{n+1}$
 and associated with homogeneous Neumann boundary conditions:
 \begin{equation}\label{bc}
\hat{u}^{'}(x)=\hat{v}^{'}(x)=0, \hskip .25cm x \in \partial \Omega.
\end{equation}
Then, we apply the Sinc-collocation method for solution
 of the system of equation (\ref{ode1})with given boundary conditions.\\
 Since the Sinc basis functions in (\ref{22}) do not have a
 derivative at endpoints $a$ and $b$, and because the Neumann boundary
 conditions (\ref{bc}) must be handled at $a$ and $b$ by approximate
 solutions,
Thus the approximate solutions for $\hat{u}(x)$ and  $\hat{v}(x)$ in Eqs(\ref{ode1}) are represented by formula \\
\begin{equation}\label{app}
\left\{\begin{array}{c} \hat{u}(x)\approx
\hat{u}_{m}(x)=c_{-N-1}w_a(x)+
\sum_{k=-N}^{N}c_{k}\xi(x)S_{k}(x)+c_{N+1}w_b(x)
,\\\\
\hat{v}(x)\approx \hat{v}_{m}(x)=\rho_{-N-1}w_a(x)+
\sum_{k=-N}^{N}\rho_{k}\xi(x)S_{k}(x)+\rho_{N+1}w_b(x),
\end{array}
                                                 \right.
\end{equation}
$$m=2N+3$$
where  $\xi(x)=(a-x)(b-x)$, we note that the first derivative of
the modified Sinc basis function of $\xi(x)S_{k}(x)$ is defined at
endpoints $a$ and $b$ and is equal to zero. Moreover, the
functions $w_a(x)$ and $w_b(x)$ satisfy in the following boundary
conditions
$$w_a(a)=w_b(b)=1,w_a^{'}(a)=w_b^{'}(b) =0,w_a(b)=w_b(a)=0, w_a^{'}(b)=w_b^{'}(a)=0,$$
and are obtained by Hermite interpolation as
$$w_a(x)=\frac{(2x+b-3a)(b-x)^2}{(b-a)^3},\hskip .5cm w_b(x)=\frac{(-2x+3b-a)(x-a)^2}{(b-a)^3}.$$

 The $4N+6$ unknown coefficients $c_{k}$ and $\rho_{k}$ in relation
(\ref{app}) are determined by substituting $\hat{u}_{m}(x)$ and
$\hat{v}_{m}(x)$ into (\ref{ode1}) and evaluating the result at
the Sinc points (\ref{24}). \\
Setting
\begin{equation}\label{app2}
 \frac{d^{i}}{d \phi ^{i}}[S_{k}(x)]=S_{k}^{(i)}(x),
\hskip 1cm 0\leq i\leq 2,
\end{equation}
 and noting that
\begin{equation}\label{app3}
\frac{d}{d x}[S_{k}(x)]=S_{k}^{(1)}(x)\phi ^{'}(x),
\end{equation}
\begin{equation}\label{app4}
\frac{d^{2}}{d x^{2}}[S_{k}(x)]=S_{k}^{(2)}(x)[\phi
^{'}(x)]^{2}+S_{k}^{(1)}(x)\phi ^{''}(x),
\end{equation}
and
\begin{equation}\label{app5}
\delta _{kj}^{(l)}= h^{l}\frac{d^{l}}{d
\phi^{l}}[S_{k}(x)]_{x=x_{j}}.
\end{equation}
By given approximations in (\ref{app}) and using
(\ref{app2}-\ref{app4}), we have
\begin{gather}
\label{appde}\hat{u}_{m}^{'}(x)=c_{-N-1}w_a^{'}(x)+
\sum_{k=-N}^{N}c_{k}\big[\xi^{'}(x)S_{k}(x)+\xi(x)\phi
^{'}(x)S_{k}^{(1)}(x)\big]+c_{N+1}w_b^{'}(x),\\
%\end{equation}
%\begin{equation}
\label{appde2}\hat{v}_{m}^{'}(x)=\rho_{-N-1}w_a^{'}(x)+
\sum_{k=-N}^{N}\rho_{k}\big[\xi^{'}(x)S_{k}(x)+\xi(x)\phi
^{'}(x)S_{k}^{(1)}(x)\big]+\rho_{N+1}w_b^{'}(x),\\
%\end{equation}
%\begin{multline}
 \nonumber \hat{u}_{m}^{''}(x)=c_{-N-1}w_a^{''}(x)+
\sum_{k=-N}^{N}c_{k}\Big[\xi^{''}(x)S_{k}(x)+2\xi^{'}(x)\phi
^{'}(x)S_{k}^{(1)}(x)+\xi(x)\phi ^{''}(x)S_{k}^{(1)}(x)\\
\label{appde1}+\xi(x)(\phi ^{'}(x))^{2}S_{k}^{(2)}(x)\Big]+c_{N+1}w_b^{''}(x),\\
%\end{multline}
%\begin{multline}\label{appde3}
\nonumber \hat{v}_{m}^{''}(x)=\rho_{-N-1}w_a^{''}(x)+
\sum_{k=-N}^{N}\rho_{k}\Big[\xi^{''}(x)S_{k}(x)+2\xi^{'}(x)\phi
^{'}(x)S_{k}^{(1)}(x)+\xi(x)\phi ^{''}(x)S_{k}^{(1)}(x)\\
\label{appde3}+\xi(x)(\phi
^{'}(x))^{2}S_{k}^{(2)}(x)\Big]+\rho_{N+1}w_b^{''}(x).
\end{gather}
 Now by substituting each terms of (\ref{ode1}) with given
approximations in (\ref{app}),(\ref{appde1}) and (\ref{appde3})
and evaluating the result at the Sinc points $x_{j}$, we can
obtain the discrete Sinc-collocation system of nonlinear equations
to  determining
 the unknown coefficients $ \{c_{k}\}_{k=-N-1}^{N+1}$ and $ \{\rho_{k}\}_{k=-N-1}^{N+1}$ as

\begin{equation}\label{diss}
\left\{\begin{array}{c} d_u \Big\{ c_{-N-1}w_a^{''}(x_j)+
\sum_{k=-N}^{N}c_{k}\Big[\xi^{''}(x_j)S_{k}(x_j)+2\xi^{'}(x_j)\phi
^{'}(x_j)S_{k}^{(1)}(x_j)+
\\\\\xi(x_j)\phi^{''}(x_j)S_{k}^{(1)}(x_j) +\xi(x_j)(\phi
^{'}(x_j))^{2}S_{k}^{(2)}(x_j)\Big]+c_{N+1}w_b^{''}(x_j)\Big\}+f\Big(c_{-N-1}w_a(x_j)+\\\\
\sum_{k=-N}^{N}c_{k}\xi(x_j)S_{k}(x_j)+c_{N+1}w_b(x_j),\rho_{-N-1}w_a(x_j)+
\sum_{k=-N}^{N}\rho_{k}\xi(x_j)S_{k}(x_j)+\\\\\rho_{N+1}w_b(x_j)\Big)-\frac{1}{\Delta
t}\Big(c_{-N-1}w_a(x_j)+
\sum_{k=-N}^{N}c_{k}\xi(x_j)S_{k}(x_j)+c_{N+1}w_b(x_j)\Big)=-\frac{u^{n}}{\Delta
t} ,
\\\\
 d_v \Big\{ \rho_{-N-1}w_a^{''}(x_j)+
\sum_{k=-N}^{N}\rho_{k}\Big[\xi^{''}(x_j)S_{k}(x_j)+2\xi^{'}(x_j)\phi
^{'}(x_j)S_{k}^{(1)}(x_j)+\\\\
\xi(x_j)\phi^{''}(x_j)S_{k}^{(1)}(x_j) +\xi(x_j)(\phi
^{'}(x_j))^{2}S_{k}^{(2)}(x_j)\Big]+\rho_{N+1}w_b^{''}(x_j)\Big\}+g\Big(c_{-N-1}w_a(x_j)+\\\\
\sum_{k=-N}^{N}c_{k}\xi(x_j)S_{k}(x_j)+c_{N+1}w_b(x_j),\rho_{-N-1}w_a(x_j)+
\sum_{k=-N}^{N}\rho_{k}\xi(x_j)S_{k}(x_j)+\\\\ \rho_{N+1}w_b(x_j)\Big)-\frac{1}{\Delta
t}\Big(\rho_{-N-1}w_a(x_j)+
\sum_{k=-N}^{N}\rho_{k}\xi(x_j)S_{k}(x_j)+\rho_{N+1}w_b(x_j)\Big)=-\frac{v^{n}}{\Delta
t},
\end{array}
                                                 \right.
\end{equation}
\hskip 3cm $j =-N-1, -N+1,...,N+1.$\\
Also, we set
$$\xi^{(i)}(x_j)=\xi^{(i)}_j,\hskip .25cm \phi^{(i)}(x_j)=\phi^{(i)}_j,
\hskip .25cm w_a^{(i)}(x_j)=w_{aj}^{(i)}, \hskip .25cm
w_b^{(i)}(x_j)=w_{bj}^{(i)},\hskip .15cm i=0,1,2 .$$ Using above
relations and (\ref{app5}), system of equations (\ref{diss}) can
be represented as
\begin{equation}\label{diss1}
\left\{\begin{array}{c} d_u \Big\{ c_{-N-1}w_{aj}^{''}+
\sum_{k=-N}^{N}c_{k}\Big[\xi^{''}_j\delta
_{kj}^{(0)}+\frac{2}{h}\xi^{'}_j\phi ^{'}_j\delta _{kj}^{(1)}+
\frac{1}{h}\xi_j\phi^{''}_j\delta _{kj}^{(1)}
+\\\\\frac{1}{h^2}\xi_j(\phi ^{'}_j)^{2}\delta
_{kj}^{(2)}\Big]+c_{N+1}w_{bj}^{''}\Big\}+f\Big(c_{-N-1}w_{aj}+
\sum_{k=-N}^{N}c_{k}\xi_j\delta
_{kj}^{(0)}+\\\\c_{N+1}w_{bj}\hskip .1cm ,\hskip .15cm
\rho_{-N-1}w_{aj}+ \sum_{k=-N}^{N}\rho_{k}\xi_j\delta
_{kj}^{(0)}+\rho_{N+1}w_{bj}\Big)-\\\\\frac{1}{\Delta
t}\Big(c_{-N-1}w_{aj}+ \sum_{k=-N}^{N}c_{k}\xi_j\delta
_{kj}^{(0)}+c_{N+1}w_{bj}\Big)=-\frac{u^{n}}{\Delta t} ,
\\\\
 d_v \Big\{ \rho_{-N-1}w_{aj}^{''}+
\sum_{k=-N}^{N}\rho_{k}\Big[\xi^{''}_j\delta
_{kj}^{(0)}+\frac{2}{h}\xi^{'}_j\phi ^{'}_j\delta _{kj}^{(1)}+
\frac{1}{h}\xi_j\phi^{''}_j\delta _{kj}^{(1)}
+\\\\\frac{1}{h^2}\xi_j(\phi ^{'}_j)^{2}\delta
_{kj}^{(2)}\Big]+\rho_{N+1}w_{bj}^{''}\Big\}+g\Big(c_{-N-1}w_{aj}+
\sum_{k=-N}^{N}c_{k}\xi_j\delta
_{kj}^{(0)}+\\\\c_{N+1}w_{bj}\hskip .1cm ,\hskip .15cm
\rho_{-N-1}w_{aj}+ \sum_{k=-N}^{N}\rho_{k}\xi_j\delta
_{kj}^{(0)}+\rho_{N+1}w_{bj}\Big)-\\\\\frac{1}{\Delta
t}\Big(\rho_{-N-1}w_{aj}+ \sum_{k=-N}^{N}\rho_{k}\xi_j\delta
_{kj}^{(0)}+\rho_{N+1}w_{bj}\Big)=-\frac{v^{n}}{\Delta t},
\end{array}
                                                 \right.
\end{equation}
\hskip 3cm $j =-N-1, -N+1,...,N+1.$\\
To obtain a
matrix representation of the equations (\ref{diss1}). We define the $m\times m$  diagonal matrix as follow:\\
\[\textbf{D}\Big(s(x)\Big)=\begin{pmatrix}
 s(x_{-N-1}) & 0 & 0 &\ldots & 0 \\
 \quad&  \quad& \quad& \quad& \quad\\
0 & s(x_{-N}) & 0 &\ldots & 0 \\
\quad&  \quad& \quad& \quad& \quad\\
\vdots& \quad&\ddots& \quad&\vdots\\
\quad&  \quad& \quad& \quad& \quad\\
0 & \ldots & 0 & \quad & s(x_{N+1})
\end {pmatrix}.\]\\
By using the above definitions  and notations in (\ref{28}), the
system (\ref{diss1})
  can be represented by  the following  form:\\
\begin{equation}\label{mat}
\left\{\begin{array}{c} d_u \big [\textbf{w}_a^{''}, \textbf{A}
,\textbf{w}_a^{''}\big ]\textbf{c}+f\Big([\textbf{w}_a, \textbf{B}
, \textbf{w}_b]\textbf{c},[\textbf{w}_a,
\textbf{B},\textbf{w}_b]\textbf{d}\Big)-\frac{1}{\Delta t}\big
[\textbf{w}_a, \textbf{B},\textbf{w}_b
\big]\textbf{c}=-\frac{\textbf{u}^{n}}{\Delta t},\\\\
d_v \big [\textbf{w}_a^{''}, \textbf{A} ,\textbf{w}_a^{''}\big
]\textbf{d}+g\Big([\textbf{w}_a, \textbf{B} ,
\textbf{w}_b]\textbf{c},[\textbf{w}_a,
\textbf{B},\textbf{w}_b]\textbf{d}\Big)-\frac{1}{\Delta t}\big
[\textbf{w}_a, \textbf{B},\textbf{w}_b
\big]\textbf{d}=-\frac{\textbf{v}^{n}}{\Delta t}
\end{array}
                                                 \right.
\end{equation}
where $\textbf{A}$ and $\textbf{B}$ are $(2N+3)\times (2N+1)$
matrices and
$\textbf{c},\textbf{d},\textbf{w}_a^{''},\textbf{w}_b^{''},\textbf{w}_a$
and $\textbf{w}_b$  are $m$-vectors  as:

$$\textbf{A}=\frac{1}{h^2}\textbf{D}\big(\xi (\phi
^{'})^{2}\big)\textbf{I}^{(2)}+\frac{1}{h}\textbf{D}\big(2\xi^{'}\phi
^{'}+\xi \phi ^{''}\big)\textbf{I}^{(1)}
+\textbf{D}(\xi^{''})\textbf{I}^{(0)},\hskip .5cm
\textbf{B}=\textbf{D}(\xi)\textbf{I}^{(0)}$$
\\
 $$\textbf{w}_a=\left(
\begin{array}{c}
 w_a(x_{-N-1}) \\
  w_a(x_{-N}) \\
  \vdots \\
w_a(x_{N+1} )
\end{array}
\right), \hskip 1cm \textbf{w}_b=\left(
\begin{array}{c}
 w_b(x_{-N-1}) \\
  w_b(x_{-N}) \\
  \vdots \\
w_b(x_{N+1} )
\end{array}%
\right), \hskip 1cm \textbf{c}= \left(
\begin{array}{c}
  c_{-N-1} \\
  c_{-N} \\
  \vdots \\
  c_{N+1} \\
\end{array}%
\right),$$ $$ \textbf{d}= \left(
\begin{array}{c}
  \rho_{-N-1} \\
  \rho_{-N} \\
  \vdots \\
  \rho_{N+1} \\
\end{array}%
\right), \hskip 1cm \textbf{w}_a^{''}=\left(
\begin{array}{c}
 w_a^{''}(x_{-N-1}) \\
  w_a^{''}(x_{-N}) \\
  \vdots \\
w_a^{''}(x_{N+1} )
\end{array}
\right), \hskip 1cm \textbf{w}_b^{''}=\left(
\begin{array}{c}
 w_b^{''}(x_{-N-1}) \\
  w_b^{''}(x_{-N}) \\
  \vdots \\
w_b^{''}(x_{N+1} )
\end{array}%
\right),$$\\
also, the matrices $\textbf{I}^{(0)}, \textbf{I}^{(1)}$ and
$\textbf{I}^{(2)}$ are $(2N+3)\times (2N+1)$ and defined in
(\ref{28}). we rewrite the system (\ref{mat}) as following form\\
\begin{equation}\label{mat1}
 \left(%
\begin{array}{cc}
   \textbf{A}_u &  \textbf{0} \\
   \textbf{0} &  \textbf{A}_v\\
\end{array}%
\right)\left(%
\begin{array}{c}
   \textbf{c} \\
   \textbf{d} \\
\end{array}%
\right)+\left(%
\begin{array}{c}
  f\big(\textbf{M}\textbf{c},\textbf{M}\textbf{d}\big) \\
 g\big(\textbf{M}\textbf{c},\textbf{M}\textbf{d}\big) \\
\end{array}%
\right)=-\frac{1}{\Delta t}\left(%
\begin{array}{c}
  \textbf{u}^{n}\\
  \textbf{v}^{n} \\
\end{array}%
\right),
\end{equation}
where
$$\textbf{A}_u=d_u \big [\textbf{w}_a^{''}, \textbf{A}
,\textbf{w}_a^{''}\big ],\hskip .5cm \textbf{A}_v=d_v \big
[\textbf{w}_a^{''}, \textbf{A} ,\textbf{w}_a^{''}\big ],\hskip
.5cm \textbf{M}=\big[\textbf{w}_a, \textbf{B} ,
\textbf{w}_b\big].$$
 The system (\ref{mat1}) is a nonlinear system of equations
which consists of $4N+6$ equations and $4N+6$ unknowns. By solving
this system by means of Newton's method , we can
 obtain  approximate solutions $u_{m}(x)$ and $v_{m}(x)$ of (\ref{ode1}) from
 (\ref{app}). Of course with initial guess zero and fix point algorithm we can obtain a starting value for Newton's method.



\section{Numerical results}\label{sec4}
In this section, by applying
the Sinc collocation method on (\ref{conEq}) to three test problems verify the
analytical results discussed in the previous sections. In all of
the examples considered in this paper, we choose $\alpha=1$ and
$d=\frac{\pi}{2}$ which yield $h=\frac{\pi}{\sqrt{2N}}$, also the
errors are reported on uniform grids
\begin{equation}\label{re}
 U=\{z_{0},z_{1},\ldots, z_{k}\},\hskip .25cm z_{s}=\frac{s}{k},\hskip .25cm s=0,1,\ldots,k.
\end{equation}
 \textbf{Example 1.} We consider a predator-prey model with Beddington-DeAngelis functional response  as  following form
\begin{equation}\label{ex1}
 \left\{{\begin{array}{l} \frac{\partial u}{\partial t}= d_u
\frac{\partial ^{2}u}{\partial x ^{2}}+u(\delta-u)-\frac{\beta
uv}{1+pu+qv},
\\\\
\frac{\partial v}{\partial t}= d_v \frac{\partial ^{2}v}{\partial
x ^{2}}+v(1-\frac{v}{u}),\hskip 2cm (x,t) \in (0,2\pi) \times
(0,\infty)
                                                 \end{array}}
                                                 \right.
\end{equation}
%\begin{equation}\label{cond2}
subject to the homogeneous Neumann boundary condition
$$\frac{\partial u}{\partial x }(0,t)=\frac{\partial u}{\partial x
}(2\pi,t)=\frac{\partial v}{\partial x
}(0,t)=\frac{\partial v}{\partial x
}(2\pi,t)=0,\hskip .75cm t \in  (0,\infty).$$
 All parameters appearing in model (\ref{ex1}) are assumed to be positive
constants. For this example, a unique positive equilibrium point
(coexistence of predator and prey) has been derrived as
$$u_s=v_s=\frac{r\delta-1-\beta+\sqrt{(1+\beta-r\delta)^2+4r\delta}}{2r}, \hskip .25cm r=p+q.$$
The qualitative analysis  for this
example is discussed in \cite{jia}. Here, we represent numerical simulations by Sinc-collocation method  for  various  values of
parameters in (\ref{ex1}). Figures \ref{pic1-1} and \ref{pic1-2} display numerical simulations of example 1 for fixed values of parameters
$N=32, p=0.6, q=0.4, \beta=12, d_v=22, \Delta t=0.01$ and for various values of $\delta$ and $d_u$. The example has been solved at final time $t=30$.
These parameters are taken from the literature \cite{jia}. In figure \ref{pic1-1} we choose $\delta=6$ such that the conditions of Theorem 1(2)(i) in \cite{jia} are satisfied. Then, by Theorem 1(2)(i) in \cite{jia}, we know that $(u_s,v_s)$ is locally asymptotically stable. In figure \ref{pic1-2} we choose $\delta=10$ such that the conditions of Theorem 1(3) in \cite{jia} are satisfied. Then, by Theorem 1(3) in \cite{jia}, we know that $(u_s,v_s)$ is unstable and non constant steady states appear.\\
The exact solution to this problem is unknown. So the accuracy of
its numerical solution will be computed using double mesh
principle, therefore for each $\delta$ the
 maximum point wise errors are estimated as
\begin{equation} \label{ex12}
 E^{m}_{i,\delta}=\max_{s}|\hat{u}_{i}^m(z_{s})-\hat{u}_{i}^{2m}(z_{s})|,
  \hskip .5 cm i=1,2,
\end{equation}
where $ \hat{u}_{1}^m=\hat{u}_m(x), \hat{u}_{2}^m=\hat{v}_m(x)$ are the
 approximation solutions of (\ref{app})  for number of  $m$ Sinc points. let
\begin{equation} \label{ex13}
 E^{m}_{\delta}=\max_{i}E^{m}_{i,\delta},
  \hskip .5 cm i=1,2.
\end{equation}

Table 1 represent values of $E^{m}_{\delta}$ for various values of $\delta$ and $N$
also for fixed values of $\beta=12, d_u=0.8, d_v=22, p=0.6, q=0.4, \Delta t=0.01$ and $ t=30$.
These results show that the errors decrease with increasing $N$.

\begin{table}[h]
\centering
  \caption{{\footnotesize Errors of $E^{m}_{\delta}$
   for example 1 with $\beta=12, d_u=0.8, d_v=22, p=0.6, q=0.4, \Delta t=0.01$ and $ t=30$}}
\begin{tabular}{|l||l||l||l||l|}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  $\delta \rightarrow $
   & $6$ & $10$  & $18$      \\

  $N \downarrow $&  &  &     \\
  \hline
  $8$ & $2.4\times 10^{-2}$ & $3.79\times 10^{-2}$ & $1.99\times 10^{-1}$  \\
  $16$ & $2.91\times 10^{-3}$ & $6.50\times 10^{-3}$ & $2.10\times 10^{-2}$   \\
  $32$ & $6.71\times 10^{-4}$ & $6.13\times 10^{-4}$ & $6.14\times 10^{-3}$   \\
  $64$ & $1.68\times 10^{-4}$ & $1.73\times 10^{-4}$ & $5.40\times 10^{-4}$  \\
  \hline
\end{tabular}
\end{table}
\begin{figure}
\begin{center}
\includegraphics[scale=.60]{figex1-1}
\\
\caption{\footnotesize{ Numerical solutions  of
example 1 for values of $\delta=6, p=0.6,q=0.4, \beta=12, d_u=0.8, d_v=22,$ with $\Delta t=0.01 $
at $t=30$ and initial conditions $(u_0,v_0)=(3,2)$ also $u_s=v_s=0.7720$.} }\label{pic1-1}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=.550]{figex1-2}
\\
\caption{\footnotesize{ Numerical solutions  of
example 1 for values of $\delta=10, p=0.6,q=0.4, \beta=12, d_u=0.8, d_v=22,$ with $\Delta t=0.01 $
at $t=30$ and initial conditions $(u_0,v_0)=(2.5+0.001\cos(\frac{x}{2}),2+0.001\cos(x))$ also $u_s=v_s=2$.} }\label{pic1-2}
\end{center}
\end{figure}
\textbf{Example 2.} Second problem  is given by
\begin{equation}\label{ex2}
 \left\{{\begin{array}{l} \frac{\partial u}{\partial t}= d_u
\frac{\partial ^{2}u}{\partial x ^{2}}+u(1-u)-\frac{
\sqrt{u}v}{1+\lambda\sqrt{u}},
\\\\
\frac{\partial v}{\partial t}= d_v \frac{\partial ^{2}v}{\partial
x ^{2}}+\beta v\Big(-\frac{\gamma+\delta v}{1+v}+\frac{
\sqrt{u}}{1+\lambda\sqrt{u}}\Big),\hskip 2cm (x,t) \in (0,\pi) \times
(0,\infty)
                                                 \end{array}}
                                                 \right.
\end{equation}
with Neumann boundary conditions and initial conditions
$$u_x(0,t)=v_x(0,t)=u_x(\pi,t)=v_x(\pi,t)=0, \hskip .5cm u(x,0)=u_0, v(x,0)=v_0. $$
The system (\ref{ex2}) is a spatial model exhibiting herd behavior in terms of square root of prey population and hyperbolic mortality $\frac{\gamma+\delta v}{1+v}$ of predator.\\
Existence and uniqueness of positive equilibrium for this model have been proved in \cite{tang}  analytically.
Here, we present some numerical simulations by Sinc procedure, also, the example has been solved at final time $t=2000, 1000$ and $t=500$.\\
Fix $\delta=0.5, \gamma=0.2, \lambda=1.5, d_u=0.01, d_v=0.8.$ According to Theorem 3.1 of \cite{tang}, the positive equilibrium $(u_s,v_s)$ is asymptotically stable for $\beta>\beta_0\approx 0.28$ and unstable for $\beta<\beta_0$. So, the system (\ref{ex2}) undergoes a Hopf bifurcation at the bifurcation value $\beta=\beta_0$. The figures \ref{pic2-1} and \ref{pic2-2}  verify the  results of Theorem 3.1 in  \cite{tang}.\\
Moreover, for the parameter values $\delta=0.5, \gamma=0.2, \lambda=1.5, d_u=0.01, d_v=5$ and $ \beta=4$, there exists a stable spatially inhomogeneous steady states. So, the system (\ref{ex2}) undergoes a pitchfork bifurcation. The figure \ref{pic2-3}  verify the  results of Theorem 3.1 in  \cite{tang}.\\
 As the exact solution $u(x,t)$ and $v(x,t)$  of (\ref{ex2}) are unknown, therefore
  the errors are estimated as (\ref{ex12}) and
  (\ref{ex13}) in example 1. Table 2
displays the results for this example for various values of $N$ and $\beta$.
\begin{table}[h]
\centering
  \caption{{\footnotesize Errors of $E^{m}_{\beta}$
   for example 2 with $\lambda=1.5, \gamma=0.2, \delta=0.5, d_u=0.01, d_v=0.8,  \Delta t=0.01$ and $ t=100$}}
\begin{tabular}{|l||l||l||l||l| }
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  $\beta \rightarrow $
   & $0.28$ & $0.32$  & $1.5$      \\

  $N \downarrow $&  &  &     \\
  \hline
  $8$ & $1.17\times 10^{-2}$ & $1.15\times 10^{-2}$ & $1.19\times 10^{-2}$  \\
  $16$ & $2.30\times 10^{-3}$ & $3.50\times 10^{-3}$ & $1.15\times 10^{-3}$   \\
  $32$ & $9.00\times 10^{-4}$ & $8.84\times 10^{-4}$ & $9.14\times 10^{-4}$   \\
  $64$ & $2.08\times 10^{-4}$ & $1.23\times 10^{-4}$ & $4.10\times 10^{-4}$  \\
  \hline
\end{tabular}
\end{table}
\begin{figure}
\begin{center}
\includegraphics[scale=.60]{figex2-1}
\\
\caption{\footnotesize{ The positive equilibrium $(u_s,v_s)=(0.4153,0.7410)$ of system (\ref{ex2}) is asymptotically stable. Here we set $\delta=0.5, \beta=0.32,\gamma=0.2, \lambda=1.5, d_u=0.01, d_v=0.8,$ and $N=32$ with $\Delta t=0.1 $
at $t=2000$ and initial conditions $(u_0,v_0)=(u_s+0.02,v_s+0.02)$.} }\label{pic2-1}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=.60]{figex2-2}
\\
\caption{\footnotesize{ Numerical simulation  of
spatially homogeneous stable periodic solution bifurcating from the unstable equilibrium $(u_s,v_s)=(0.4153,0.7410)$ of system (\ref{ex2}). Here we set $\delta=0.5, \beta=0.28,\gamma=0.2, \lambda=1.5, d_u=0.01, d_v=0.8,$ and $N=32$ with $\Delta t=0.1 $
at $t=1000$ and initial conditions $(u_0,v_0)=(u_s+0.02,v_s+0.02)$. } }\label{pic2-2}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=.60]{figex2-3}
\\
\caption{\footnotesize{ Numerical simulation  of pitchfork bifurcation of positive constant equilibrium of system (\ref{ex2}) for parameter values as $\delta=0.5, \beta=4,\gamma=0.2, \lambda=1.5, d_u=0.01, d_v=5 ,$ and $N=32$ with $\Delta t=0.01 $
at $t=500$ and initial conditions $(u_0,v_0)=(u_s+0.2\cos(x),v_s+0.05\cos(x))$. The positive constant equilibrium $(u_s,v_s)=(0.4153,0.7410)$ is unstable.
 } }\label{pic2-3}
\end{center}
\end{figure}


\textbf{Example 3.} We consider prey-predator reaction-diffusion system \cite{dim}
\begin{equation}\label{ex3}
 \left\{{\begin{array}{l} \frac{\partial u}{\partial t}= d_u
\frac{\partial ^{2}u}{\partial x ^{2}}+R u(1-\frac{u}{K})-\beta \frac{
u^2 v}{1+\lambda u^2},
\\\\
\frac{\partial v}{\partial t}= d_v \frac{\partial ^{2}v}{\partial
x ^{2}}+\gamma \frac{
u^2 v}{1+\lambda u^2}-\delta v,\hskip 2cm (x,t) \in (-20,20) \times
(0,\infty)
                                                 \end{array}}
                                                 \right.
\end{equation}
subject to boundary  conditions
$$u_x(-20,t)=v_x(-20,t)=u_x(20,t)=v_x(20,t)=0, $$
and initial conditions
$$ u(x,0)=1-\frac{1}{2}\sin^{10}(\pi(x-20)/40),\hskip .5cm v(x,0)=\frac{1}{4}\sin^{10}(\pi(x-20)/40).$$
The growth rate of the prey $f(u)=R u(1-\frac{u}{K})$ is logistic and the predator's  functional response
$\frac{u^2 }{1+\lambda u^2}$ is Holling type III,
the ratio $\gamma/\beta$ and parameter $R$ represent the maximal \emph{per capita} predator and prey birth rates respectively, $\delta$ is the \emph{per capita} predator death rate and $K$ is the prey carrying capacity. \\ In the following, we adopt the parameters of system (\ref{ex3}) according to the corresponding parameters in \cite{dim},  our simulations confirm the results obtained in \cite{dim}.
Figures present simulations of prey and predator for various cases of these parameters. In these figures the diffusion coefficient associated with the prey is much greater
than the diffusion coefficient of the predators($d_u=1, d_v=10^{-5}$). Figure \ref{pic3-1}(a) illustrates the prey-predator interaction characterized by the type III
(sigmoid) functional response, in which the rate of attack of the predator ($v$) accelerates at first and then decelerates towards satiation. Such sigmoid functional responses are typical of natural enemies which readily switch from one
food species to another and/or which concentrate their feeding in areas where certain resources are most abundant.
These are called prey-dependent responses because the feeding rate of consumers is dependent only on the density of
prey.\\
 This characteristic profile for the prey population density appears to be affected when we increase its carrying
capacity $K$ and the maximal \emph{per capita} prey birth rate $R$ (see figure \ref{pic3-1}(b), or increase the prey carrying capacity  $K$ together with decreasing the\emph{ per capita} predator death rate $\delta$ (see figure \ref{pic3-2}). Also, according to (\ref{ex12}) and (\ref{ex13}), values of error for this example are given in table 3 for various values of $N$ and $K$.
\begin{table}[h]
\centering
  \caption{{\footnotesize Values of error
  for example 3 for various values of $N$ and $K$ with $\lambda=10, \gamma=0.001, \delta=0.05, \beta=1, R=0.075, d_u=1, d_v=10^{-5}, \Delta t=0.01$ and $t=100$.}}
\begin{tabular}{|l||l||l||l||l| }
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  $K \rightarrow $
   & $1$ & $5$  & $10$      \\

  $N \downarrow $&  &  &     \\
  \hline
  $8$ & $5.21\times 10^{-1}$ & $3.61\times 10^{-1}$ & $8.27\times 10^{-2}$  \\
  $16$ & $4.18\times 10^{-2}$ & $2.18\times 10^{-2}$ & $9.01\times 10^{-3}$   \\
  $32$ & $6.51\times 10^{-3}$ & $5.01\times 10^{-3}$ & $7.81\times 10^{-4}$   \\
  $64$ & $8.01\times 10^{-4}$ & $6.12\times 10^{-4}$ & $3.13\times 10^{-4}$  \\
  \hline
\end{tabular}
\end{table}
\begin{figure}
\begin{center}
\includegraphics[scale=.7]{figex3-1}
\\
\caption{\footnotesize{ Comparative results on the densities of the prey $u(x,t)$ and predators $v(x,t)$ of example 3 for values of $d_u=1, d_v=10^{-5}, \delta=0.01, \beta=100,\gamma=0.1, \lambda=10,$ and $N=32$ with $\Delta t=0.01 $
at $t=200$.} }\label{pic3-1}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[scale=.65]{figex3-2}
\\
\caption{\footnotesize{ Comparative results on the densities of the prey $u(x,t)$ and predators $v(x,t)$ of example 3 for values of $d_u=1, d_v=10^{-5}, R=0.075, \beta=1,\gamma=0.001, \lambda=10,$ and $N=32$ with $\Delta t=0.01 $
at $t=200$.} }\label{pic3-2}
\end{center}
\end{figure}

\section{Conclusions}
In this article, a numerical method was employed successfully for
solving reaction-diffusion systems arising from two-species ecological models.
This approach  is based on the implicit Euler method for
temporal discretization and the Sinc collocation method in the
spatial direction. Results from numerical experiments indicate the efficiency
and accuracy of proposed method. Also, from the figures
\ref{pic1-1}-\ref{pic3-2}  of numerical simulations, we get some useful
information about the biological behaviors of species. Comparisons in given
examples show the agreement of the approximate solutions with those presented in  \cite{jia,tang,dim}.

%\section{Bibliography}

\begin{thebibliography}{99}
\bibitem{can}  Cantrell R., Cosner C.,{\em Spatial ecology via reaction-diffusion equations}, Wiley series in mathematical
 and computational biology. Wiley ,Chichester (2003)
\bibitem{hal}  Holmes E., Lewis M., Banks J.,  Veit R., {\em Partial differential equations in ecology: spatial interactions
and population dynamics}, Ecology, 75(1), 17–29, (1994).
\bibitem{Tur}  Turing A. M., {\em The chemical basis of morphogenesis, Phil}. Trans.R. Soc.B,  237, 37-72, (1952),
\bibitem{van} vande Koppel J.,  Crain C.M.,{\em Scale-dependent inhibition drives regular tussock spacing in a fresh water marsh},Am. Nat. ,168, 136-147,(2006).
\bibitem{duf} Dufiet V., Boissonade J., {\em Dynamics of Turing pattern monolayers close to onset},Phys. Rev. , 53,4883-4892,(1996).
\bibitem{sun}Sun G.-Q., Zhang G., Jin Z., Li L.,{\em Predator cannibalism can give rise to regular spatial pattern in a predator-prey system}, Nonlinear Dynam. , 58
75-84,(2009).
\bibitem{sun1}Sun G.Q., Jin Z., Li L., Li B.L., {\em Self-organized wave pattern in a predator-prey model}, Nonlinear Dynam, , 60, 265-275,(2010).

\bibitem {liu} Liu Q.X.,  Sun G.Q., Jin Z., Li B.L.,{\em Emergence of spatiotemporal chaos arising from far-field break up of spiral waves in the plankton ecological systems},
Chin.Phys.B,,18, 506-515,(2009).
\bibitem {Wan} Wang W., Cai Y., Wu M., Wang K., Li Z. , {\em Complex dynamics of a reaction-diffusion epidemic model}, Nonlinear Anal. RWA, ,13, 2240-2258,(2012).
\bibitem{baek}  Baek H. ,  Jung DI.,  Wang Z., {\em Pattern formation in a semi-ratio-dependent predator-prey system with diffusion}, Discr
Dyn Natur Soc.  doi:10.1155/2013/657286,(2013).
\bibitem{garv}  Garvie M, {\em Finite-difference schemes for reaction-diffusion equations modeling predator-pray interactions in
MATLAB}. Bullet MathBiol.,69, 931-56,( 2007).
\bibitem{Mi}  Mimura M. , Sakaguchi H.,  Matsushita M.,{\em Reaction-diffusion modelling of bacterial colony patterns}, Physica A. ,
282,283-303,(2000).
\bibitem{pet}  Petrovskii S. , Morozov AY.,  Venturino  E., {\em Allee effect makes possible patchy invasion in a predator-prey system}, Ecol
Lett , 5, 345-52,(2002).
\bibitem{jia} Jiang H. ,  Wang L., Yao R., {\em Numerical simulation and qualitative analysis for apredator-prey
model with B-D functional response}, Mathematics and Computers in
Simulation , 117, 39-53,(2015)
\bibitem{wan} Wang  W.,  Liu QX.,  Jin Z., Spatiotemporal complexity of a ratio-dependent predator-prey system. Phys Rev E. ,
75, 1539-3755,(2007).
\bibitem{garv1} Garvie M.,  Trenchea C., {\em Spatiotemporal dynamics of two generic predator-prey models}. J Biol Dyn.,4, 559-570,(2010).
\bibitem{kol}  Kolade M. Owolabi,  Kailash C. Patidar, {\em Higher-order time-stepping methods for time-dependent
reaction-diffusion equations arising in biology}, Applied
Mathematics and Computation, 240, 30-50,( 2014)
 \bibitem{fat}Fatemeh Shakeri, Mehdi Dehghan, {\em The finite volume spectral element method to solve Turing models in the
biological pattern formation}, Computers and Mathematics with
Applications,  62, 4322-4336,(2011)
\bibitem{mar} Marcus R. Garvie , Philip K. Maini , Catalin Trenchea, {\em An efficient and robust numerical algorithm for estimating parameters
in Turing systems}, Journal of Computational Physics
, 229 (19),7058-7071,(2010).
\bibitem{kol1} Kolade M.Owolabi
and Kailash C. Patidar, {\em Numerical simulations of multicomponent
ecological models with adaptive methods}, Theoretical Biology and
Medical Modeling , DOI10.1186/s12976-016-0027-4, (2016)
\bibitem{mar1} Marcus R. Garvie , John Burkardt , Jeff Morgan,
{\em Simple Finite Element Methods for Approximating Predator-Prey
Dynamics in Two Dimensions Using Matlab}, Bull Math
Biol, 77,548-578,(2015)
\bibitem{st} F. Stenger, {\em A Sinc Galerkin method of solution of boundary value
problems}, Math. Comp. , 33,  35-109,(1979).
\bibitem{lu}  Lund J.,    Bowers K., {\em Sinc Methods for Quadrature and Differential Equations}, SIAM, Philadelphia,
PA,(1992).
\bibitem{st1}  Stenger F., {\em Numerical Methods Based on Sinc and Analytic Functions}. Springer, New York,(1993).
\bibitem{bi} Bialecki B. , {\em Sinc-collocation methods for two-point
boundary value problems}. IMA J. Numer. Anal.,11, 357-375,(1991).
\bibitem{RMA} Rashidinia J. ,  Nabati M.,  Barati A., {\em Sinc-Galerkin method for solving nonlinear weakly
singular two point boundary value problems}, International Journal
of Computer Mathematics, , DOI: 10.1080/00207160.2015.1085027,(2015)
\bibitem{rash}  Rashidinia J.  ,   Barati A.,  Nabati M., {\em Application of Sinc-Galerkin method to singularly
perturbed parabolic convection-diffusion problems}, Numer. Algor. , 66,643-662,(2014).
\bibitem{ga} M. EL-GAMEL,  {\em The sinc-Galerkin method for solving singularly-perturbed reaction-diffusion problems}, Electronic Transactions on Numerical
Analysis. , 23,  129-140,(2006).
\bibitem{to} Okayama T., Matsuo T.,
Sugihara  M., {\em Sinc-collocation methods for weakly singular
Fredholm integral equations of the second kind}, Journal of
Computational and Applied Mathematics., 234, 1211-1227,(2010).
\bibitem{sa}  Saadatmandi A., Dehghan M., {\em The use of Sinc-collocation method for solving multi-point
boundary value problems}, Commun Nonlinear Sci Numer Simulat ,
17(2), 593-601,(2012).
\bibitem{tang}  Tanga X., Song Y.,{\em Bifurcation analysis and Turing instability in a diffusive predator-prey model with herd behavior and hyperbolic mortality}
,Chaos, Solitons and Fractals , 81, 303-314,(2015)

\bibitem{dim}  Narcisa A.,   Dimitriu G., {\em On a prey-predator reaction-diffusion system with Holling type III functional response}, Journal of computational and applied mathematics , 235(2), 366-379,(2010).
\end{thebibliography}

\end{document}