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

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

\begin{document}
\title[\hfilneg EJDE-2018/183\hfil 
 Elliptic sectors and Euler discretization]
{Elliptic sectors and Euler discretization}

\author[N. Mohdeb, A. Fruchard, N. Mehidi \hfil EJDE-2018/183\hfilneg]
{Nadia Mohdeb, Augustin Fruchard, Noureddine Mehidi}

\address{Nadia Mohdeb \newline
Laboratoire de Math\'ematiques Appliqu\'ees,
Universit\'e A. Mira, Bejaia, Alg\'erie}
\email{n\_mohdeb@hotmail.com}

\address{Augustin Fruchard \newline
 Laboratoire de Math\'ematiques, Informatique et  Applications,
Universit\'e de Haute Alsace, Mulhouse, France}
\email{augustin.fruchard@uha.fr}

\address{Noureddine Mehidi \newline
Laboratoire de Math\'ematiques Appliqu\'ees\\
Universit\'e A. Mira, Bejaia, Alg\'erie}
\email{manogha@yahoo.fr}

\dedicatory{Communicated by Adrian Constantin}

\thanks{Submitted September 1, 2018. Published November 14, 2018.}
\subjclass[2010]{34C25, 34C37, 34A34, 39A05}
\keywords{Elliptic sector; nonhyperbolic equilibrium point;
 homoclinic orbit;
\hfill\break\indent  S-invertible; Euler discretization}

\begin{abstract}
 In this work we are interested in the elliptic sector of
 autonomous differential systems with a degenerate equilibrium point
 at the origin, and in their Euler discretization.
 When the linear part of the vector field at the origin has two zero
 eigenvalues, then the differential system has an elliptic sector,
 under some conditions.  We describe this elliptic sector and we show
 that the associated Euler discretized system has an elliptic sector
 converging to that of the continuous system when the step size of the
 discretization tends to zero.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\allowdisplaybreaks

\section{Introduction}

In this work we consider the planar differential system
\begin{equation}
\begin{gathered}
\dot{x}=  P(x,y) \\
\dot{y}=  Q(x,y)
\end{gathered}  \label{1}
\end{equation}
where $P$ and $Q$ are analytic functions from $\mathbb{R}^2$ to
$\mathbb{R}$. Also we consider their Euler discretization
\begin{equation}
\begin{gathered}
x_{n+1}=x_n+hP(x_n,y_n) \\
y_{n+1}=y_n+hQ(x_n,y_n)
\end{gathered}   \label{discretise}
\end{equation}
where $h>0$ is the step size of the discretization.
We assume that the origin is an isolated equilibrium point of \eqref{1}
and  that system \eqref{1} has an elliptic sector $S_0$.

The main aim of this work is to explore to what extent the discrete
system  \eqref{discretise} presents also an elliptic sector $S_{h}$
and whether $S_{h}$ tends to $S_0$ in the sense of the Hausdorff
distance, when $h$ tends to zero. \newline
In order to define an elliptic sector of \eqref{1}, we  consider a
circle $C$ with center  $(0,0)$ and radius $r$, containing no other
equilibria than the origin. We will assume that there exist
two solutions $\gamma _1$
and $\gamma _2$ of system \eqref{1} tending to the origin;
we assume for example that $\gamma _1$ tends to the origin when  $t$
tends to $+\infty $ and  $\gamma _2$ tends to the same point when $t$
tends to $-\infty $.
We denote by $\gamma _1^{\ast }$ and $\gamma _2^{\ast }$ the
respective corresponding orbits to the solutions
$\gamma _1^{\ast }$ and $\gamma _2^{\ast }$.
Let $M_1$ and $M_2$ be the respective intersection points of
$\gamma _1^{\ast }$ and $\gamma _2^{\ast }$ with the  circle
$C$ such that, taking into account the direction of the orbits,
$M_1$ is the last intersection point of $\gamma _1^{\ast }$
with the circle $C$ and $M_2$ is the first intersection  point of
$\gamma_2^{\ast }$ with $C$. Let $\sigma $ be the closed curve made
up of the two segments $OM_1$ and $OM_2$ parts of the two orbits
$\gamma_1^{\ast }$ and $\gamma _2^{\ast }$,  the origin and  the
oriented arc of $C$ joining  $M_1$ to $M_2$ in the forward direction.
The region $R_{\sigma }$ delimited  by $\sigma $ is said to be a
\textit{sector}. An \textit{elliptic sector} is a sector containing only
nested homoclinic or parts of nested homoclinic orbits
(Figure \ref{sectellip}) \cite{Andronov}.
We recall that a solution is said to be homoclinic if it is defined on
$\mathbb{R}$  entirely and tends to the origin when $t$
tends to $+\infty $ as well as when $t$ tends to $-\infty $.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=5cm, height=3.2cm]{fig1a} % secteur.png
\quad
\includegraphics[width=4.cm, height=3.3cm]{fig1b} % elliptique.png
\end{center}
\caption{A sector and an elliptic  sector}
\label{sectellip}
\end{figure}

To pose the problem raised above correctly, we must define the notion
of homoclinic orbits of system \eqref{discretise}, even when the map
$$
F_{h}{:}( x,y) \mapsto (x+hP(x,y),y+hQ(x,y))
$$
is not invertible. This leads us to define the  S-invertibility notion.
The map $F_{h}$ is  said to be S-invertible if, for any $m_0$ in
$\mathbb{R}^2$, there exist two reals $R>0$ and $h_0>0$ such that
for any $h$ in $]0,h_0]$, there exists a unique $m_{-1,h}$ in $\mathbb{R}^2$
for which $\| m_{-1,h}\| \leq R$ and $F_{h}(m_{-1,h})=m_0$.
In this case,  we have necessarily
\[
\underset{h\to 0}{\lim }~m_{-1,h}=m_0.
\]
In other words, $F_{h}$ is S-invertible if any point of the  plane has a unique
``good predecessor'' by the application $F_{h}$. The notation S comes from
 nonstandard analysis (NSA): in the  context of NSA, a planar map is
said to be  S-invertible if any limited point of the plane has  a
unique limited predecessor; for the standard functions, this notion
coincides with the usual invertibility.

Let $U$ be an open set of the disc $D(0; r_0)$ and  simply connected.
We say that $U$ is an elliptic sector of the discrete system
\eqref{discretise} if any solution emanating from $U$ is homoclinic.

The works of  Beyn \cite{bey87 le bon}. Fiedler and  Scheurle
\cite{Fied96} and that of  Zou \cite{Zou2003} deal with the problem of
the persistence of the non-degenerate homoclinic orbits of an autonomous
system after the discretization of the latter. They give an error estimate
of order $O(h^{d})$, for the difference between the homoclinic solution
of the differential equation and that of the associated discrete equation,
where $h$ is the step size of the method of discretization and $d$ its order.
They also give the length $l(h)$ of the parameter interval over which the
homoclinic orbit persists. On the other hand, given an autonomous
differential system which has an equilibrium  point at $(0,0)$ not
necessarily hyperbolic and an unstable center manifold $W_c^{u}$,
it is shown in \cite{beyn87ab} that, in a small neighborhood of the origin,
the discrete  system associated by a one step method has, under some
conditions, an invariant manifold close to $W_c^{u}$.
In the hyperbolic case, it is shown in  \cite{beyn87} that, under some
conditions, the phase  portrait of the differential system is correctly
reproduced in the associated discretization by a one step method,
on an arbitrary time  interval. In \cite{Eirola3}, the author
examines the conditions which make the solution of the differential system 
close to that of the associated discretized one on an infinite time  interval. 
A similar study has been made in \cite{Garay1} for the structurally stable 
systems without periodic solutions. On the other hand, it has been 
shown in \cite{Feckan2} that in the hyperbolic case, the local stable
 manifold of the  discrete system tends to that of the corresponding 
differential system when the step size tends to zero 
(see also \cite{Feckan4,Garay5}). An extension of this result to nonautonomous 
differential systems is given in \cite{Garay6}. A Taylor expansion  
approximation of this manifold is given in \cite{Eirola1}.
In \cite{Feckan3} it is shown that, in the hyperbolic case,
the maps defining the vector fields and the associated discretized
system by the one step method are uniformly topologically equivalent.
Other results about numerical calculus of homoclinic, heteroclinic
and periodic orbits are established in  
\cite{beyn87b}-\cite{beyn2002}, \cite{Doed89}-\cite{Fried93},
\cite{Zou97} and \cite{Zou2003}.

It is known that, when system \eqref{1} has an elliptic  sector,
then the associated Jacobian matrix  $M$ of the function
$(x,y)\mapsto (P(x,y),Q(x,y))$ at point $(0,0)$, has two zero eigenvalues
(cf. \cite[p. 241]{lefschetz}); the origin is a non-hyperbolic equilibrium
point. Two situations are possible. The first one is that where the
matrix $M$ is null. In this case, the behavior of the solutions near
the origin is very  complex;
If the smallest degree of the nonlinear terms of \eqref{1} is $m$,
then the neighborhood of the origin will be split into $2(m+1)$
parabolic, hyperbolic or elliptic sectors. The number of elliptic
sectors existing in system \eqref{1} depends on the index
 of the equilibrium point (cf. \cite{perko}, p 151): let $C$ be a
Jordan curve containing $(0,0)$ and no other critical point of \eqref{1}
in its interior; then the index of the equilibrium point $(0,0)$ with
respect to \eqref{1} is given by
\[
I_{\eqref{1}}(0,0) = I_{\eqref{1}}(C)
=\frac{1}{2\pi}\oint_C \frac{PdQ-QdP}{P^2+Q^2}.
\]
The second situation is that where the matrix $M$ can  be reduced by a
linear transformation to  $M'=\begin{pmatrix}
0 & 1 \\
0 & 0
\end{pmatrix}$; it is in this situation that we are interested in this work.
In this  case,  system \eqref{1} is reduced by a sequence of analytic
changes of variables to the form (see \cite{Andronov})
\begin{equation}
\begin{gathered}
\dot{x}=  y \\
\dot{y}=  ax^{r}\big(1+k(x)\big)+bx^{p}y\big(1+g( x) \big)
+y^2f( x,y)
\end{gathered}   \label{norm}
\end{equation}
where $f$, $g$ and $k$ are analytic functions, such that
$k(0)=g(0)=f( 0,0) =0$, $a$ and $b$ are real  parameters,
$r$ and $p$ are integer parameters, satisfying  $r=2m+1$,
$m\geq 1$, $a<0$, $b\neq 0$, $p\geq 1$, $p$ odd and either

\begin{itemize}
\item $p=m$ and $\Delta =b^2+4(m+1)a\geq 0$, or

\item  $p<m$.
\end{itemize}

Without loss of generality, by a  linear change of variable,
 we assume in all the following that $ b $ is strictly positive.
In all this work, we place ourselves in the situation $p=m=1$
and $\Delta =b^2+8a>0$.  These assumptions are relatively natural.
Indeed, the assumption $p=m=1$ is  generic in  some sense.
The assumption $\Delta\geq0$ is necessary to get an elliptic sector,
and the assumption $\Delta\neq0$ is generic.
The affine  transformation
$(x,y)\mapsto (x,\sqrt{-a}y)$ together with the time change
$\tau =\sqrt{-a}t$ and the parameter change $b\mapsto b/\sqrt{-a}$,
allow us to assume that $a=-1$. Thus we  will be interested in the
elliptic sector of system
\begin{equation}
\begin{gathered}
\dot{x}=  y \\
\dot{y}=  -x^{3}( 1+k( x) ) +bxy\big( 1+g(x) \big) +y^2f( x,y)
\end{gathered}   \label{sytcontmodif}
\end{equation}
and in its persistence in the associated Euler discretized system
\begin{equation}
\begin{gathered}
x_{n+1}=x_n+hy_n \\
y_{n+1}=y_n+h\Big(-x_n^{3}\big(1+k(x_n)\big)+bx_ny_n\big(1+g(x_n)
\big)+y_n^2f( x_n,y_n) \Big)
\end{gathered}  \label{sytcont33}
\end{equation}
The main result of this work is the following.

\begin{theorem} \label{theorx3copy1}
There exists $h_0>0$ such that for any $h$ in $]0,h_0]$, there exists
a subregion $S_{h}$ of the elliptic sector $S_0$ of
system \eqref{sytcontmodif} with the following properties:
\begin{itemize}
\item Any solution of system \eqref{sytcont33} starting from
 $S_{h}$ is homoclinic.

\item When $h$ tends to zero, $S_{h}$
tends to $S_0$ in the sense of the Hausdorff distance.
\end{itemize}
\end{theorem}

We do not try to obtain a maximal region of homoclinic orbits
of \eqref{sytcont33}, whose structure may be highly complicated.
In particular we do not study the behavior of orbits of \eqref{sytcont33}
starting close to the boundary of $S_0$.
A more precise study of these orbits remains to be done.
The speed of convergence of $S_h$ to $S_0$ as $h$ tends to $0$ is another
interesting problem which is not discussed here.

The article  is organized as follows.
In section two, we give a local  description of the elliptic sector
of system \eqref{sytcontmodif} whose existence is stated above.
To do so, we will use the  results obtained during the study of the
global behavior of the solutions of the model example
\begin{equation}
\begin{gathered}
\dot{x}=  y \\
\dot{y}=  -x^{3}+bxy
\end{gathered}   \label{sytcontx3}
\end{equation}
In section three, we give a complete proof of theorem \ref{theorx3copy1}.
In the last section we deal with the elliptic sector of the family
of differential systems which are diffeomorphic to system
\eqref{sytcontmodif} and with its persistence in the corresponding
discretized systems.

\section{Description of the continuous system}

We are first interested in the global behavior of the solutions of system
\eqref{sytcontx3} (Figure \ref{1Champs}).
The following proposition describes the phase portrait of \eqref{sytcontx3}.
The symmetry with respect to the $y$-axis allows to consider only the
solutions of \eqref{sytcontx3} starting at points with positive abscissas.
Let us denote $\Delta=b^2-8$ and $\alpha _1=( b-\sqrt{\Delta }) /4$
and $\alpha _2=( b+\sqrt{\Delta }) /4$ the solutions of the equation
$-2\alpha ^2+b\alpha -1=0$.

\begin{proposition}\label{propx3}
Let $( x_0,y_0) $ be a point of the plane such that
$x_0>0$ and let $\gamma(t,x_0,y_0) $ be the solution  of system
\eqref{sytcontx3} emanating from this point.
\begin{itemize}
\item[(1)] If $y_0<\alpha _1x_0^2$, then the solution
 $\gamma(t,x_0,y_0) $ is homoclinic.
\item[(2)] If $\alpha _1x_0^2<y_0<\alpha _2x_0^2$, then
\begin{gather*}
\lim_{t\to -\infty}x(t)=\lim_{t\to -\infty}y(t)=0,\\
\lim_{t\to +\infty}x(t)=\lim_{t\to +\infty }y(t)=+\infty
\end{gather*}
Moreover,
\[
\lim_{t\to -\infty}\frac{y(t)}{x(t)^2}=\alpha _1
\quad\text{and}\quad \lim_{t\to +\infty}\frac{y(t)}{x(t)^2}
=\alpha _2.
\]

\item[(3)] If $\alpha _1x_0^2<y_0$, then the trajectory $\gamma $
crosses the $y$-axis at a point whose y-coordinate is positive. Also
\[
\lim_{t\to +\infty}x(t)=\underset{t\to \pm
\infty }{\lim }y(t)=-\lim_{t\to -\infty}x(t)=+\infty .
\]
\end{itemize}
\end{proposition}

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=6cm, height=5cm]{fig2} % phasexy.png
\end{center}
\caption{The phase portrait of system \eqref{sytcontx3}.}
\label{1Champs}
\end{figure}

Hereafter, we give only a sketch of the proof since it is similar to proof
of the proposition \ref{propecl}.
Note however  that it uses the traditional tools of trap  trajectories,
taking into account the fact that the curves $y=\alpha _1x^2$
and $y=\alpha _2x^2$ are solutions of system \eqref{sytcontx3},
the blow-up theory \cite{Pelletier} and other tools of qualitative
study of differential equations.  Indeed,  the  application
for system \eqref{sytcontx3} of the two quasi-homogeneous blow-ups
of degree $(1,2)$
\begin{equation}
x=u,\quad y=u^2v\quad \text{with }d\tau =udt,  \label{ecl11}
\end{equation}
and
\begin{equation}
x=wz, \quad y=z^2\quad \text{with }d\eta =zdt  \label{ecl12}
\end{equation}
produce respectively the systems
\begin{equation}
\begin{gathered}
u'=  uv \\
v'=  -1+bv-2v^2
\end{gathered}   \label{syschgt}
\end{equation}
and
\begin{equation}
\begin{gathered}
w'=  1+{\dfrac{1}{2}}w^{4}-\dfrac{b}{2}w^2 \\
z'=  -\dfrac{1}{2}w^{3}z+\dfrac{b}{2}wz
\end{gathered}   \label{systecl}
\end{equation}
The phase portraits of systems \eqref{syschgt} and \eqref{systecl}
are illustrated in Figure \ref{eclatement1}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm, height=6cm]{fig3a} % eclat1.png
\includegraphics[width=6cm, height=6cm]{fig3b} % eclat2.png
\end{center}
\caption{Phase portraits of systems \eqref{syschgt} and
\eqref{systecl} for $b=3$.}
\label{eclatement1}
\end{figure}

The interpretation of these phase portraits allows us to prove the
stated result.

Now, we prove the following proposition giving a local description
of the elliptic  sector of system \eqref{sytcontmodif}.

\begin{proposition}\label{propecl}
System \eqref{sytcontmodif} has an elliptic sector. Moreover,
the behavior of the solutions of \eqref{sytcontmodif} near the
origin is similar to that of system \eqref{sytcontx3}.
More precisely, in a neighborhood of the origin, solutions of
\eqref{sytcontmodif} starting below the parabola
$y=\alpha _1x^2$ are homoclinic and solutions starting above
the parabola $y=\alpha _2x^2$ are not homoclinic.
\end{proposition}

\begin{proof}
The two previous blow-ups applied to \eqref{sytcontmodif}
give respectively the two systems
\begin{equation}
\begin{gathered}
u'=  uv \\
v'=  -1+bv-2v^2-k(u)+bvg(u)+uv^2f( u,u^2v)
\end{gathered}  \label{syteclat6}
\end{equation}
and
\begin{equation}
\begin{gathered}
w'=  1+\frac{w^{4}}{2}(1+k(wz))-\frac{bw^2}{2}(1+g(wz))
-\frac{wz}{2}f(wz,z^2) \\
z'=  -\frac{w^{3}z}{2}(1+k(wz))+\frac{bwz}{2}(1+g(wz))
+\frac{z^2}{2}f(wz,z^2)
\end{gathered}  \label{syteclat2}
\end{equation}
The blow-up $x=XY$, $y=-Y^2$ together with $d\xi =Ydt$ applied to
system \eqref{sytcontmodif} give
\begin{equation}
\begin{gathered}
X'=  -1-\frac{X^{4}}{2}(1+k(XY))-\frac{bX^2}{2}(1+g(XY))
+\frac{XY}{2}f(XY,-Y^2) \\
Y'=  \frac{X^{3}Y}{2}(1+k(XY))+\frac{bXY}{2}(1+g(XY))
-\frac{Y^2}{2}f(XY,-Y^2)
\end{gathered}   \label{syteclat3}
\end{equation}
System \eqref{syteclat6} has an unstable node on $A_1=(0,\alpha _1)
$ and a saddle on $A_2=(0,\alpha _2)$ (Fig. \ref{portrait1}).
The $v$-axis is invariant under  system \eqref{syteclat6}.
The continuity of the solutions with respect to the initial
conditions gives the existence of $\delta >0$ such that for any
$u_0$ in the interval $]-\delta ,\delta ] $,
the solution of \eqref{syteclat6} starting at $(u_0,0)$ tends to
$A_2$ when $\tau$ tends to $-\infty $ and enters in the lower
half-plane after some time $\tau_0 >0$ (depending on $u_0$).

\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm, height=6cm]{fig4a} % portr.png
\includegraphics[width=6cm, height=6cm]{fig4b} % eclat2p.png
\end{center}
\caption{the phase portraits of system \eqref{syteclat6} near
 $A_1$ and $A_2$ and of system \eqref{syteclat2}
near $B_1$, $B_2$, $B_3$ and $B_4$ for $f(x,y)=g(x)=k(x)=x$ and $b=3$.}
\label{portrait1}
\end{figure}

System \eqref{syteclat2} has as equilibria (Figure \ref{portrait1})
two saddles on
$B_1=(-\sqrt{2\alpha _1},0)$ and on
$B_2=(\sqrt{2\alpha _1},0)$ and two nodes, stable on
$B_3=(-\sqrt{2\alpha _2},0)$ and unstable on
$B_4=(\sqrt{2\alpha _2},0)$. The solutions in the neighborhoods of
$B_1$, $B_2$, $B_3$ and
$B_4$ having the tangents $w=\sqrt{2\alpha _1}$, $w=-\sqrt{2\alpha _1}$,
$w=\sqrt{2\alpha _2}$ or $w=-\sqrt{2\alpha _2}$, correspond in the
$( x,y) $ plane to the solutions of system \eqref{sytcontmodif}
whose  order two approximations near $(0,0)$ are given by
$y=\alpha_1x^2$ and $y=\alpha _2x^2$. \newline
The $w$-axis is invariant under system \eqref{syteclat2}.
By continuity of the solutions of  \eqref{syteclat2} with respect to the
initial conditions, those that lie near $B_1$ and $B_2$
sufficiently close to the lines $w=-\sqrt{2\alpha _2}$\ and
$w=\sqrt{2\alpha _2}$, with the initial conditions $(w_0,z_0)$
such that $-\sqrt{2\alpha _2}<w_0<\sqrt{2\alpha _2}$,
move close to the $w$-axis, then cross the $z$-axis on points whose
$y$-coordinates are close to zero.

In the plane $(x,y)$, this means that the corresponding solutions lie
in the neighborhood of the origin close to the curves $y=\alpha _1x^2$
and $y=\alpha _2x^2$, then cross the $y$-axis near the origin.\newline
System \eqref{syteclat3} has no equilibrium point and the $X$-axis
is invariant under this system. When $Y=0$,  system \eqref{syteclat3} becomes
\begin{equation}
\begin{gathered}
X'= -1-\dfrac{1}{2}X^{4}-\dfrac{b}{2}X^2 \\
Y'=  0
\end{gathered}   \label{sysecl44}
\end{equation}
By continuity of solutions of \eqref{syteclat3} with respect to the initial
conditions, there exist some  solutions emanating from  points close enough
to the $X$-axis, that cross the $Y$-axis on points whose y-coordinates are
sufficiently close to zero.

We deduce that, in the lower half-plane, there is no solution of system
\eqref{sytcontmodif} which  tends to $(0,0)$.
A simultaneous interpretation of results of the three blow-ups
applied to system \eqref{sytcontmodif} allows us to conclude.
\end{proof}

\section{Homoclinic orbits of the discrete system \eqref{sytcont33}}

In this part we give the proof of theorem \ref{theorx3copy1}.
This proof is based on the properties  of solutions  of system
\eqref{sytcontx3}.

For convenience, we have chosen to use the formalism of Nonstandard Analysis
(NSA for short) \cite{rob} in the IST version due to E. Nelson \cite{nelson};
In particular, we use the tools of stroboscopy \cite{Sari}
and of permanence principle \cite{Isambert}.
In a classical framework, several statements can play the role of the
permanence principle of NSA, e.g.\ the Kaplun's extension theorem \cite{kaplun}
(see also \cite[page 27]{lagerstrom})
or  \cite[Theorem 2.2.2]{eckhaus}.

From now on, the functions $f,g,k$ of system
\eqref{sytcontx3} are assumed to be standard, and $h>0$ is fixed infinitesimal.
By the Transfer Principle of NSA, it is sufficient to prove that,
for all limited $(x_0,y_0)\in S_0$ in the $S$-interior of $S_0$,
the solution $( x_n,y_n)_{n\in \mathbb{Z}}$
of \eqref{sytcont33} starting from  $(x_0,y_0)$ is homoclinic.
We recall that $(x_0,y_0)$ is in the $S$-interior of $S_0$ if there exists a
standard $r>0$ such that the disk of center $(x_0,y_0)$ and of radius $r$
is in $S_0$.
For $n<0$, the points $( x_n,y_n)$ are given
recursively by Lemma \ref{lemx3fgh4}.

The proof of Theorem \ref{theorx3copy1} essentially uses the idea that
for $\alpha _2\lnsim a\lnsim \alpha _1$, the region
$\{(x,y)\in \mathbb{R}^2;y<ax^2\}$ is positively invariant under
system \eqref{sytcontx3}. By
$a\lnsim b$ we mean that $a$ is less than $b$ and not infinitly close to $b$.

To prove Theorem \ref{theorx3copy1}, we need a second lemma which
is inspired by the idea above. The proof of this lemma requires to
highlight an intermediate discrete system. Put
$\ell ( x,y)=-x^{3}(1+k(x))+bxy(1+g(x))+y^2f(x,y)$ and
$F_{h}( x,y) =(x+hy,y+h\ell (x)) $.

Also there exists a function $\psi$ defined near $x=0$, satisfying
 $\ell(x,\psi(x))=0$ and
\[
\psi (x)=\frac{1}{b}x^2+o(x^2)\text{.}
\]
Then, the component $\dot{y}$ of the vector field \eqref{sytcontmodif}
 becomes zero on the curve $y=\psi(x)$.

Now, let $(x_0,y_0)$ be a limited point in the elliptic sector of
system \eqref{sytcontmodif}, not infinitely close to its border.
Several cases are possible, but  only one is presented here:
$x_0>0$ and $0<y_0<\psi (x_0)$. The other cases are similar.

The change of variables $x=\varepsilon X$, $y=\varepsilon ^2Y$ and
the time change $\tau =\varepsilon t$ where $\varepsilon \simeq 0$,
transform  system \eqref{sytcontmodif} into
\begin{equation}
\begin{gathered}
X'=  Y \\
Y'=  -X^{3}(1+k(\varepsilon X))+bXY(1+g(\varepsilon
X))+\varepsilon Y^2f(\varepsilon X,\varepsilon ^2Y)
\end{gathered}   \label{sytcontloupfgh}
\end{equation}
where $(')=d/d\tau $. By using the short shadow lemma
\cite{Diener}, it appears that the solutions of system
\eqref{sytcontloupfgh} starting at limited points $(X_0,Y_0)$
have for shadows in a limited time interval,  solutions of system
\eqref{sytcontx3}. Let $( \bar{x}(t),\bar{y}(t)) $ denote the solution
of system \eqref{sytcontmodif} starting at point $(x_0,y_0)$, and let
$(x_n,y_n)_{n\in \mathbb{N}}$ be the solution of system
\eqref{sytcont33} started from the same point.

For $nh$ bounded, by using the stroboscopy lemma \cite{Sari},
\begin{equation}
(x_n,y_n)\simeq (\bar{x}(nh),\bar{y}(nh))  \label{stroboscop}
\end{equation}
By permanence \cite{Diener}, there exists $N$ in $\mathbb{N}$ such that
$Nh$ is infinitely large and  for any $n<N$, approximation \eqref{stroboscop}
remains true. Then
\[
(\bar{x}(Nh),\bar{y}(Nh))\simeq (0,0)
\]
Let $\varepsilon =-\bar{x}(Nh)$ and put $x_n=\varepsilon X_n$,
$y_n=\varepsilon ^2Y_n$ and $\bar{h}=h\varepsilon $. System
\eqref{sytcont33} becomes
\begin{equation}
\begin{gathered}
X_{n+1}=X_n+\bar{h}Y_n \\
\begin{aligned}
Y_{n+1}&=Y_n+\bar{h}\Big(-X_n^{3}(1+k(\varepsilon
X_n))+bX_nY_n(1+g(\varepsilon X_n)) \\
&\quad +\varepsilon Y_n^2f(\varepsilon X_n,\varepsilon ^2Y_n)\Big)
\end{aligned}
\end{gathered}\label{sysdiscloupe11}
\end{equation}

The following lemmas are used in the proof of theorem \ref{theorx3copy1}.

\begin{lemma} \label{lemx3fgh3}
 Let $\alpha_2\lnsim a \lnsim \alpha_1$. Any solution
$( X_n,Y_n) _{n\in \mathbb{N}}$ of system \eqref{sysdiscloupe11}
starting at a limited point $( X_0,Y_0) $ in the region
$\{(x,y)\in \mathbb{R}^2/x<0,\;y< a x^2\}$ remains in it.
\end{lemma}

\begin{proof}
We use induction on $n$:
\[
X_{n+1}=YX_n+h_n \leq X_n+ha(X_n)^2=X_n(1+haX_n)
\]
As $X_n$ is assumed negative, then $X_{n+1}<0$.
Otherwise,
\begin{align*}
Y_{n+1}-aX_{n+1}^2
&=( Y_n-aX_n^2) ( 1+\bar{h}( b-2a+bg(\varepsilon X_n)) X_n) \\
&\quad +  \bar{h}( -2a^2+ba-1+bag(\varepsilon X_n)-k(\varepsilon
X_n)) X_n^{3} \\
&\quad +\bar{h}\varepsilon Y_n^2f(\varepsilon X_n,
\varepsilon ^2Y_n)-a \bar{h}^2Y_n^2
\end{align*}
On the other hand, $-2a^2+ba-1\gnsim0$. Since $g$ and $k$ are standard,
continuous and  satisfy $g(0)=k(0)=0$, we deduce that
\[
-2a^2+ba-1+bag(\varepsilon X_n)-k(\varepsilon X_n)\gnsim 0.
\]
Hence
\begin{equation}
\Big(-2a^2+ba-1+bag(\varepsilon X_n)-k(\varepsilon X_n)\Big)
X_n^{3}+\varepsilon Y_n^2f(\varepsilon X_n,\varepsilon ^2Y_n)<0
\label{terme}
\end{equation}
Thus $Y_{n+1}-aX_{n+1}^2<0$, when $Y_n-aX_n^2<0$.
\end{proof}

\begin{lemma} \label{lemx3fgh4}
Any limited point in the plane has a unique limited
predecessor by system \eqref{sytcont33}.
\end{lemma}

\begin{proof}
With the previous notation $\ell$ and $F_h$, we must show that for any
limited point  $(u,v)$ of the plane, there exists a unique limited point
$(x,y)$ such that
\begin{equation}
F_{h}(x,y)=(u,v).  \label{Fh}
\end{equation}
We mention  that such a point, if it exists, is necessarily infinitely
close to $(x,y)$. Moreover, $x$ is given by $x=u-hy$.
Then, it suffices to show that there exists a unique $y\in ] v-1,v+1[$
satisfying
\begin{equation}
v=y+hl(u-hy,v).  \label{pred11}
\end{equation}
For a fixed $( u,v) $, we denote
\[
L( y) =v-y-hl(u-hy,v)
\]
Since $l(u-h(v+1),v)$ is limited,
\[
L( v+1) =-1-hl(u-h(v+1),v)<0
\]
In the same way,
\[
L( v-1) =1-hl(u-h(v-1),v)>0
\]
Also, for any limited $y$,
\[
L'( y) =-1-h( -h\frac{\partial l}{\partial x}(
u-hy,v) +\frac{\partial l}{\partial y}( u-hy,v) ) <0
\]
Since $L$ is continuous, the equation $L( x) =0$ has a unique limited
solution $y$ in $] v-1,v+1[ $.
\end{proof}

We fix a limited point $(x_0,y_0)$.
Let $( x_{-n},y_{-n}) _{n\in \mathbb{N}^\ast }$ be the predecessors
sequence of system \eqref{sytcont33} defined by lemma \ref{lemx3fgh4}.
Then, this sequence is uniquely defined, as long as $( x_{-n},y_{-n})$
is limited, by
\begin{gather*}
x_{-n}=x_{-n-1}+hy_{-n-1} \\
y_{-n}=y_{-n-1}+h\ell(x_{-n-1},y_{-n-1})
\end{gather*}
satisfying $( x_{-1},y_{-1}) \simeq ( x_0,y_0) $
and $( x_{-n-1},y_{-n-1}) \simeq ( x_{-n},y_{-n})$.

\begin{lemma} \label{lemx3fgh5}
Let $a\in \mathbb{R}$ such that $\alpha_2\lnsim a\lnsim
\alpha_1$ and $(X_0,Y_0)$ a limited point in the plane such that $
0<Y_0<aX_0^2$ and $X_0>0$. The solution $( X_{-n},Y_{-n})
_{n\in \mathbb{N}^\ast }$ of system \eqref{sysdiscloupe11} starting at $
(X_0,Y_0)$, does not leave the region of the plane
\[
\{(X,Y)\in \mathbb{R}^2:  X>0,\; 0<Y<aX^2\}
\]
\end{lemma}

\begin{proof}
We use induction. We will show that $Y_{-1}<aX_{-1}^2$;
The same reasoning allows us to show that $Y_{-n}<aX_{-n}^2$ for any $n>1$.
We denote $\delta =-2a^2+ba-1+bag(\varepsilon X_{-1})-k(\varepsilon X_{-1})$
and $\beta = \varepsilon f(\varepsilon X_{-1},\varepsilon ^2Y_{-1})-a\bar{h}$.
 We have
\begin{eqnarray*} Y_0-aX_0^2
=(Y_{-1}-aX_{-1}^2) \Big( 1+\bar{h}\big( b-2a+bg(\varepsilon X_{-1})\big)
X_{-1}\Big) + \bar{h} \delta X_{-1}^{3}+\bar{h}\beta
\end{eqnarray*}
We will distinguish  two cases:
If $\beta\geq 0$,
since $\delta>0$ we have
\[ \label{ineg1}
0\geq Y_0-aX_0^2\geq ( Y_{-1}-aX_{-1}^2)
(1+\bar{h}(b-2a+bg(\varepsilon X_{-1}))X_{-1}).
\]
This means that $Y_{-1}-aX_{-1}^2<0$.

On the other hand,
\[
\delta X_{-1}^{3}+\beta Y_{-1}^2
= \beta  Y_{-1}(Y_{-1}-aX_{-1}^2)+ a\beta X_{-1}^2\Big(Y_{-1}
 +  \frac{\delta }{a\beta}X_{-1}\Big).
\]
Since $\delta\not\simeq 0$, $\delta /(a\beta)$ is infinitely large,
then if $\beta <0$,
\[
\frac{\delta }{a\beta}X_{-1}+aX_{-1}^2
= \Big(\frac{\delta }{a\beta}+aX_{-1}\Big)X_{-1} <0.
\]
Thus,
\[
\delta X_{-1}^{3}+\beta Y_{-1}^2\geq\beta (Y_{-1}-aX_{-1}^2)(
Y_{-1}+aX_{-1}^2),
\]
hence
\begin{align*}
&Y_0-aX_0^2 \\
&\geq  ( Y_{-1}-aX_{-1}^2) \Big( 1+\bar{h}(
b-2a+bg(\varepsilon X_{-1}) ) X_{-1} + \beta(
Y_{-1}+aX_{-1}^2)\Big)
\end{align*} %\label{ineg2}
It results in this case that  $Y_{-1}<aX_{-1}^2$. \
\end{proof}

\begin{proof}[Proof of theorem \ref{theorx3copy1}]
By applying the stroboscopy lemma, the solution
$( x_n,y_n) _{n\in \mathbb{N}}$ reaches a point of the region
\[
\{(x,y)\in \mathbb{R}^2:  x<0, \;0<y<\psi (x)\}
\]
and as long as $y_n<\alpha_2x_n^2$ then $x_{n+1}<0$.
Furthermore, this solution reaches a point
$(x_{q},y_{q})$ infinitely close to the origin, in the region
\[
\{(x,y)\in \mathbb{R}^2: x<0, \psi (x)<y< a x^2\}\text{.}
\]
On the other hand, by lemma \ref{lemx3fgh3},
\[
    \forall n>q, x_n<0, \psi (x_n)<y_n< a x_n^2
\]
and the sequence $( x_n,y_n)_{n\in \mathbb{N}}$ is convergent;
its limit is the stationary point $(0,0) $ of \eqref{sytcont33}.

Likewise, by using  lemmas \ref{lemx3fgh4} and \ref{lemx3fgh5},
it results that the sequence $( x_{-n},y_{-n}) _{n\in \mathbb{N}}$
 converges to $( 0,0)$
when $n$ tends to $+\infty $.
\end{proof}

\section{Discretization of the family of planar differential systems
 diffeomorphic to \eqref{sytcontmodif}}

In this section we consider a system of the form \eqref{1} which is
diffeomorphic to system  \eqref{sytcontmodif} in the following sense:
there exists a diffeomorphism $\varphi{:} \mathbb{R}^2\to\mathbb{R}^2$
given by
\[
 \varphi =(\varphi _1,\varphi _2):(x,y)\mapsto
(u,v)=(\varphi _1(x,y),\varphi _2(x,y))
\]
such that a function from $\mathbb{R}$ to $\mathbb{R}^2$,
$t\mapsto(x(t),y(t))$ is a solution of \eqref{1} if and only if
the function  $t\mapsto(u(t),v(t))=\varphi(x(t),y(t))$ is a solution
of  \eqref{sytcontmodif} (with the  variables denoted $u$ and $v$
instead of $x$ and $y$).

It will also be  possible to make a time change, but the question which
interests us concerns the elliptic sectors and these do not depend on
the parametrization of the integral curves. We then assume immediately
that the two systems are diffeomorphic without time change.
We mention that it is not the same to discretize a differential system
before or after applying diffeomorphisms.

System \eqref{1} can be written as
\begin{equation}
\begin{gathered}
\dot{x}=  ax+by+P_1(x,y) \\
\dot{y}=  cx+dy+Q_1(x,y)
\end{gathered}\label{sytgen}
\end{equation}
where $a$, $b$, $c$ and $d$ are reals and $P_1$ and $Q_1$ are analytic
functions whose Taylor expansion near the origin begins with terms of
total degree on $x$ and $y$  greater than or equal to two.

We know that when system \eqref{sytgen} has an elliptic sector,
then the two eigenvalues of the associated matrix of its linear part
are zero  \cite{lefschetz,perko}. This means that
\begin{gather*}
a=-d, \\
a^2+bc=0 \\
|a|+|b|+|c|\neq 0
\end{gather*}
The most generic case is when $a\neq 0$, $b\neq 0$
and $c\neq 0$. In this case, the change of variables
\[
 \chi {:}(u,v) \mapsto ( x,y) =\chi ( u,v)
\]
defined by
\begin{gather*}
x=  v \\
y=  cu-av+\bar{P}(v,cu-av)
\end{gather*}
transforms system \eqref{sytgen} into system \eqref{sytcontmodif}.
The Jacobian of $\chi $ at point $(0,0)$ is ($-c$)$\neq 0$.
This means that in some neighborhood of $(0,0)$, the map $\chi $
is bijective. It  suffices to put $\varphi =\chi ^{-1}$. We write
\begin{gather*}
\varphi _1(x,y)=  \bar{f}(x,y) \\
\varphi _2(x,y)=  x
\end{gather*}
where $\bar{f}$ is an analytic function such that $\bar{f}(0,0)=0$.
It follows that the diffeomorphism $\varphi$ transforms system
\eqref{discretise} into
\begin{equation}
\begin{gathered}
u_{n+1}=u_n+hv_n \\
v_{n+1}=v_n+h l(u_n,v_n)+h^2G(x_n,y_n),
\end{gathered}  \label{sysdiscgen}
\end{equation}
where
\begin{align*}
G(u_n,v_n) &= \frac{1}{d_n}\Big(\frac{( u_{n+1}-u_n) ^2}{2}
 \frac{\partial ^2\varphi _1}{\partial u^2}(u_n,v_n)
 +\frac{( v_{n+1}-v_n) ^2}{2}
 \frac{\partial ^2\varphi _1}{\partial v^2}(u_n,v_n) \\
&\quad + (u_{n+1}-u_n)( v_{n+1}-v_n)
 \frac{\partial ^2\varphi_1}{\partial u\partial v}\Big)
 \frac{\partial \varphi _2}{\partial u}(u_n,v_n) \\
&\quad + o\Big( \| ( u_{n+1},v_{n+1}) -(u_n,v_n) \| ^{3}\Big)
\end{align*}
with
\[
 d_n=\frac{\partial \varphi _1}{
\partial u}(u_n,v_n)\frac{\partial \varphi _2}{\partial v}
(u_n,v_n)-\frac{\partial \varphi _2}{\partial u}(u_n,v_n)\frac{
\partial \varphi _1}{\partial v}(u_n,v_n)
\]
We can check that
$G$ can be written for any $(x,y)\in \mathbb{R}^2$ as
\[
G(x,y)= xyG_1(x,y)+y^2G_2(x,y)
\]
where $G_1$ and $G_2$  are analytic functions in $\mathbb{R}^2$.

In the same way as system \eqref{discretise},  system \eqref{sytcontmodif}
has an elliptic sector. Hence, system \eqref{1} has also an elliptic sector.
We will show by passing through system \eqref{sysdiscgen} obtained
by diffeomorphism from \eqref{sytcontmodif}
and by adapting the same proof as for theorem \ref{theorx3copy1},
that system \eqref{discretise} has an elliptic sector tending to that
of system \eqref{1} when $h$ tends to zero.


We remark that the study of elliptic sectors in system
(\ref{sytcontmodif}) and in its discretized one $(\ref{sytcont33})$
is similarly generalized to   system $(\ref{norm})$ and to its discretized
one obtained by Euler in the case where $p=m$ and $\Delta >0$.
The case $\Delta =0$ is treated separately.


\begin{thebibliography}{00}

\bibitem{Andronov} A. A. Andronov, E. A. Leontovich, I. I. Gordon, A. G. Maier;
\emph{Qualitative theory of second order dynamic systems},
John Wiley and Sons, New York.Toronto, 1973.

\bibitem{bey87 le bon} W.-J. Beyn;
\emph{The effect of discretization on homoclinic orbits},
in: T.\ K\"{u}pper, et al. (Eds.), Bifurcation,
Analysis, Algorithms, Applications, Birkh\"{a}user, Basel  (1987), 1-8.

\bibitem{beyn87ab} W.-J. Beyn, J. Lorenz;
\emph{Center manifolds of dynamical systems under discretization},
Numer. Funct. Anal. and Optimiz., 9(3\&4) (1987), 381-414.

\bibitem{beyn87} W.-J. Beyn;
\emph{On the numerical approximation on
phase portraits near stationary points}, SIAM J. Numer Anal, Vol 24, no. 5, 987.

\bibitem{beyn87b} W.-J. Beyn;
\emph{On invariant closed curves for one-step methods},
 Numer. Math, 51  (1987), 103-122.

\bibitem{beyn90} W.-J. Beyn;
\emph{The numerical computation of connecting orbits in dynamical systems},
IMA J. Num. Anal., 9  (1990), 379-405.

\bibitem{beyn94} W.-J. Beyn;
\emph{Numerical analysis of homoclinic orbits
emanating from a Takens-Bogdanov point}, IMA J. Num. Anal., 14 (1994), 381-410.

\bibitem{beyn97} W.-J. Beyn and J.-M. Kleinhauf (1997), The numerical
computation of homoclinic orbits of maps, SIAM J. Num. Anal. 34, 1207-1236.

\bibitem{beyn99} W.-J. Beyn, M. Stiefenhofer;
\emph{A direct approach to homoclinic orbits in the fast dynamics of
singularly perturbed systems}, J. Dyn. Differential Eq., Vol 11, no. 4  (1999),
671-709.

\bibitem{beyn2000} W.-J. Beyn, J. Schropp;
\emph{Runge-Kutta discretizations of singularly perturbed gradient equations},
BIT, Vol 40, no. 3 (2000),  415-433.

\bibitem{beyn2002} W.-J. Beyn, Barnabas M. Garay;
\emph{Estimates of variable stepsize Runge-Kutta methods for sectorial
evolution equations with nonsmooth data}, Appl. Num. Math. 41 (2002), 369-400.

\bibitem{Diener} F. Diener, G. Reeb;
\emph{Analyse non standard}, Hermann, 1989.

\bibitem{Diener2} F. Diener, M. Diener;
\emph{Nonstandard analysis in practice}, Springer Universitext, 1995.

\bibitem{Doed89a}  E. J. Doedel, M. J. Friedman;
\emph{Numerical computation of heteroclinic orbitS},
J. Comput. Appl. Math., 26  (1989),  159-170.

\bibitem{Doed89} E. J. Doedel;
\emph{Numerical computation of heteroclinic orbits},
J. Comput. Appl.  Math. 26,  (1989), 155-170.

\bibitem{Doed94} E. J. Doedel, M. J. Friedman, A. C. Monteiro;
\emph{On locating connecting orbits}, Applied Math. and Comput., 239
 (1994), 231-239.

\bibitem{Doed97} E. J. Doedel, M. J. Friedman, B. I. Kunin;
\emph{Successive continuation for locating connecting orbits}, Numerical
algorithms, 14  (1997), 103-124.

\bibitem{eckhaus} W. Eckhaus;
\emph{Asymptotic Analysis of Singular Perturbation},
Studies in Mathematics and its Applications, North Holland, 1979.

\bibitem{Eirola3} T. Eirola;
\emph{Aspects of backward error analysis of numerical ODE's},
J. of Comp. and Applied Mathematics, 45 (1993), 65-73.

\bibitem{Eirola1} T. Eirola;
\emph{Numerical Taylor expansions for invariant manifolds},
Helsinki Universiry of Technology Institute of Mathematics Research Reports,
2003.

\bibitem{Feckan2} M. Feckan;
\emph{Asymptotic behavior of stable manifolds},
Proceedings of the American Mathematical Society, Vol. 111 no.  2, 1991.

\bibitem{Feckan3} M. Feckan;
\emph{Discretization in the method of averaging}, Proceedings of the
American Mathematical Society, Vol. 113 no.  4 (1991), 1105-1113.

\bibitem{Feckan4} M. Feckan;
\emph{The relation between a flow and its discretization},
 Math. Slovaca 42 (1992), 123-127.

\bibitem{Fied96} B. Fiedler, J. Scheurle;
\emph{ Discretization of homoclinic orbits, rapid forcing and
``invisible'' chaos}, Memoirs of the AMS ,  119 (570) 79, 1996.

\bibitem{Fried91} M. J. Friedman, E. J. Doedel;
\emph{Numerical computation and continuation of invariant manifolds
connecting fixed points}, SIAM J. Num. Anal., 28 (3) (1991), 789-808.

\bibitem{Fried93} M. J. Friedman, E. J. Doedel;
\emph{Computational methods for global analysis of homoclinic and
heteroclinic orbits: A case study}, J. Dyn. and Differential Eq.,
 Vol. 5 (1993), 231-239.

\bibitem{Garay6} B. M. Garay;
\emph{Discretisation of semilinear differential equations with an
exponential dichotomy}, Computers Math. Appl., Vol. 28 no. 1-3
 (1994), 23-35.

\bibitem{Garay5} B. M. Garay;
\emph{Discretisation and some qualitative properties of ordinary
 differential equations about equilibria}, Acta Math. Univ. Comenianae,
vol. LXIII  (1994),  249-275.

\bibitem{Garay1} B. M. Garay;
\emph{Various closeness concepts in numerical ODE's},
Computers Math. Appl. Vol. 31 no. 4/5  (1996), 113-119.

\bibitem{Isambert} \'E. Isambert, V. Gautheron;
\emph{lire l'Analyse Non Standard}, Bull. Soc. Math. Belgique,
suppl\'ement ``Nonstandard Analysis''  (1996), 29-49.

\bibitem{kaplun} S. Kaplun;
\emph{Fluid Mechanics and Singular Perturbations},
P. A. Lagerstrom, L. N. Howards, C. S. Liu Ed., Academic Press, New York, 1967.

\bibitem{lagerstrom} P. A. Lagerstrom;
\emph{Matched Asymptotic Expansions: Ideas and technics},
Applied Mathematical Sciences 76, Springer Verlag, 1988.

\bibitem{lefschetz} S. Lefschetz;
\emph{Differential equations: geometric theory},
Interscience publishers, 1957.

\bibitem{nelson} E. Nelson;
\emph{Internal Set Theory: a new approach to nonstandard analysis},
Bull. Amer. Math. Soc., 83   (1977), 1165--1198.

\bibitem{Pelletier} M. Pelletier;
\emph{\'Eclatements quasi homogènes}, Annales de la facult\'e des sciences de
Toulouse $6^{e}$ s\'erie, tome 4, No. 4  (1995), 879-937.

\bibitem{perko} L. Perko;
\emph{Differential equations and dynamical systems},
Springer-Verlag, 2000.

\bibitem{rob} A. Robinson;
\emph{Nonstandard Analysis}, North-Holland, 1967.

\bibitem{Sari} T. Sari;
\emph{Stroboscopy and averaging},
Colloque trajectorien \'a la m\'emoire de G. Reeb et J. L. Callot,
Strasbourg-Obernai, 1995.

\bibitem{Zou97} Y.-K. Zou, W.-J. Beyn;
\emph{Invariant manifolds for nonautonomous systems with application
to one-step methods}, J. Dyn. and Differential Eq. Vol 10  (1997),  379-407.

\bibitem{Zou2003} Y.-K. Zou, W.-J. Beyn;
\emph{On manifolds of connecting orbits in discretizations of dynamical
 systems}, Nonlinear Analysis, 52  (2003), 1499-1520.

\end{thebibliography}

\end{document}