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

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

\begin{document}
\title[\hfilneg EJDE-2017/280\hfil Exponential stability]
{Exponential stability of solutions of nonlinear fractionally perturbed
ordinary differential equations}

\author[E. Brestovansk\'a, M. Medve\v{d} \hfil EJDE-2017/280\hfilneg]
{Eva Brestovansk\'a, Milan Medve\v{d}}

\address{Eva Brestovansk\'a \newline
Department of  Economics and Finance,
Faculty of Management,
 Comenius University, Odboj\'arov 10,
831 04 Bratislava, Slovakia}
\email{Eva.Brestovanska@fm.uniba.sk}

\address{Milan Medve\v{d} (corresponding author)\newline
Department of Mathematical Analysis and Numerical Mathematics,
Faculty of Mathematics, Physics and Informatics,
Comenius University, 842 48 Bratislava, Slovakia}
\email{Milan.Medved@fmph.uniba.sk}

\dedicatory{Communicated by Mokhtar Kirane}

\thanks{Submitted May 5, 2017. Published November 10, 2017.}
\subjclass[2010]{34A08, 34A12, 34D20, 37C75}
\keywords{Riemann-Liouville integral; Riemann-Liouville derivative;
 \hfill\break\indent Caputo derivative; fractional differential equation;
 exponential stability}

\begin{abstract}
The main aim of this paper is to prove a theorem on the exponential stability
of the zero solution of a class of integro-differential equations,
 whose right-hand sides involve the Riemann-Liouville fractional integrals
of different orders and we assume that they are polynomially bounded.
Equations of that type can be obtained e.g. from fractionally damped
pendulum equations, where the fractional damping terms depend on the
Caputo fractional derivatives of solutions. The set of initial values
of solutions that converge to the origin is also determined.
We also prove an existence and uniqueness theorem for this type of equations,
which we use in the proof of the stability theorem.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks

\section{Introduction}\label{s:1}

Recently, fractional differential equations with fractional derivatives
and fractional integrals of different types have attracted many scientists
from various disciplines  due to their wide applications.
The most known are the Riemann-Liouville and the Caputo derivatives and
differential equations with these derivatives. The basic theory of
fractional differential equations and many references can be found in
the monographs  \cite{LAK,POD,ZHO}. These derivatives are defined as follows:

The Riemann-Liouville fractional derivative of a function
$u\colon[0, \infty)\to \mathbb{R}$ of an order $\alpha \in (0, 1)$ is
\begin{equation}\label{RLD}
{}^{RL}D^\alpha u(t) := \frac{1}{\Gamma(\alpha)}\frac{d}{dt}
\int_0^t(t-s)^{\alpha - 1}u(s)ds
\end{equation}
and the Caputo derivative is
\begin{equation}\label{CD}
{}^CD^\alpha u(t) := \frac{1}{\Gamma(1- \alpha)}\int_0^t(t-s)^{-\alpha}u'(s)ds,
\end{equation}
where $u'(t) = \frac{du(t)}{dt}, \Gamma(z)
= \int_0^\infty \tau^{z -1}e^{-\tau}d\tau$ is the Euler Gamma function
and the integral
\begin{equation}\label{RLI}
^{RL}I^\alpha u(t) =  \frac{1}{\Gamma(\alpha)}\int_0^t(t-s)^{\alpha - 1}u(s)ds
\end{equation}
is called the  Riemann-Liouville fractional integral of order $\alpha$
of the function $u(t)$,
provided that the right-hand sides of \eqref{RLD}, \eqref{CD} and \eqref{RLI},
 respectively, exist.

If $u\colon [0, \infty)\to \mathbb{R}^N$,
$u(t) = \big (u_1(t), u_2(t),\dots, u_N(t)\big )$, then, we define the
Riemann-Liouville fractional derivative of the mapping $u(t)$ of order
$\alpha$ as
\begin{equation}
^{RL}D^\alpha u(t) = \big (^{RL}D^\alpha u_1(t), ^{RL}D^\alpha u_2(t),\dots,
 ^{RL}D^\alpha u_N(t)\big)
\end{equation}
and  the Riemann-Liouville fractional integral as
\begin{equation}
^{RL}I^\alpha u(t) = \big (^{RL}I^\alpha u_1(t), ^{RL}I^\alpha u_2(t),\dots,
^{RL}I^\alpha u_N(t)\big ).
\end{equation}
An influence of viscous fluids on vibrating systems is often modeled by
using the Riemann-Liouville or Caputo fractional derivative.
These derivatives play the role of damping force, called the fractional damping.
The well known Bargley-Torvik equation (see \cite{BT})
\begin{equation}
    u''(t) + A\,^CD^{\frac{3}{2}}u(t) = au(t) + \phi(t),
\end{equation}
modelling the motion of a rigid plate immersing in a viscous liquid,
is one of the equations descibing the motion with the fractional damping
term $A^CD^{\frac{3}{2}}u(t)$.

It is well known that the system of linear fractional differential equations
\begin{equation}
    D^\alpha x(t) = Ax(t),\quad x(t) \in \mathbb{R}^N,\ \alpha \in (0, 1),
\end{equation}
where $D^\alpha x(t)$ is the Riemann-Liouville or the Caputo derivative
of $x(t)$ of the order $\alpha \in (0, 1)$ and $A$ is a constant matrix,
do not have exponentially stable solutions, but asymptoticall stable only.
The equilibrium $x=0$ of this equation is
 asymptotically stable if and only if  $|\arg (\lambda)| > \frac{\alpha \pi}{2}$
for all eigenvalues of the matrix $A$. In this case all components of $x(t)$
decay towards 0  like $t^{-\alpha}$
 (see \cite{MAT1,MAT2,PET1,PET2,LZ}).

 We will show that fractional integro-differential equations of the form
\begin{equation}\label{E4}
    \dot x(t) =Ax(t) + f\big (t, x(t), ^{RL}I^{\alpha_1}x(t), \dots,
 ^{RL}I^{\alpha_m}x(t)\big), \quad x(t) \in  \mathbb{R}^N
\end{equation}
can have exponentially stable solutions, where the solution is defined
as in the next definition.

\begin{definition} \rm
A mapping $x\colon [0, T) \to \mathbb{R}^N, $ where $0 < T \leq \infty$,
is a solution of the equation \eqref{E4} satisfying the initial condition
$x(0) = x_0 \in \mathbb{R}^N$ if it is continuously differentiable on the
 interval $(0, T)$, continuous on  $[0, T)$ and it satisfies the equality
\eqref{E4} for all $t \in (0, T)$.
If $T = \infty$, then this solution is called global.
\end{definition}

Equations of the form \eqref{E4} can be obtained from the following linear
multi-fractional pendulum equation
\begin{equation}\label{I6}
u''(t) + \lambda_1(t)^CD^{\beta_1}u(t) +\dots + \lambda_m(t)^CD^{\beta_m} u(t)
 + \lambda u'(t) + \omega^2u(t) = 0
\end{equation}
with $m$ fractional and one ordinary damping terms, which can be written as
 a system of the form \eqref{E4} with
\begin{gather*}
    A= \begin{pmatrix}
    0 & 1\\
    -\omega^2 & - \lambda\end{pmatrix},\quad
    x(t)=\begin{pmatrix}
    x_1(t)\\x_2(t)\end{pmatrix}, \\
\begin{aligned}
 &f(t, x(t), ^{RL}I^{\alpha_1}x(t),\dots, ^{RL}I^{\alpha_m}x(t) )\\
&=\begin{pmatrix}
    0\\-\lambda_1(t)^{RL}I^{\alpha_1}x_2(t) -\dots
- \lambda_m(t)^{RL}I^{\alpha_m}x_2(t) - \lambda(t)x_2(t)\end{pmatrix},
\end{aligned}
\end{gather*}
where $x_1(t) = u(t)$, $x_2(t) = u'(t)$, $\beta_i \in (0, 1)$,
$\alpha_i = 1 - \beta_i$, $\lambda_i(t)$ $i=1,2,\dots, m$ are continuous
 functions on $[0, \infty)$ and $\lambda > 0$,
$\omega > 0$ are constants.

Let us recall an analysis of the fractional vibration equation
\begin{equation}\label{FVE}
    u''(t) + b^CD^\alpha u(t) + cu(t) = 0,\quad \alpha \in (0, 1),
\end{equation}
where $b > 0$, $c > 0$ are constants, given in the papers \cite{LD} and \cite{NA}.
It is proven in \cite{NA} that if $u(t)$ is a solution of the equation
\eqref{FVE} satisfying the initial conditions $u(0) = u_0$, $u'(0) = u_1$,
then its Laplace transform is
\[
\mathcal{L}[u(t)] = U(s) = \frac{s + bs^{\alpha -1} + c}{s^2 + bs^\alpha + c}u_0
+ \frac{1}{s^2 + bs^\alpha + c}u_1,
\]
the characteristic equation
\begin{equation}\label{CHE}
s^2 + bs^\alpha + c = 0
\end{equation}
 has a couple of complex conjugate roots
\[
 s_{1, 2} = \beta \pm i\sigma = r^{\pm i\Theta},\quad \beta < 0,
\quad \sigma > 0,\quad  r = \sqrt{\beta^2 + \sigma^2} > 0,\quad
\frac{\pi}{2} < \Theta < \pi
\]
and the fundamental solution $\phi_1(t)$ with
\[
\mathcal{L}[\phi_1(t)] = \Phi_1(s)
= \frac{s + bs^{\alpha -1} + c}{s^2 + bs^\alpha + c}
\]
has the form
\[
\phi_1(t) = Ce^{\beta t}\cos{\sigma t} + D e^{\beta t}\sin{\sigma t}
+  \int_0^\infty K_\alpha(\tau)e^{-\tau t}d\tau,
\]
where $C, D$ are functions of the variables $r, \Theta$. The function
\[
 f_1(t) = Ce^{\beta t}\cos{\sigma t} + D e^{\beta t}\sin{\sigma t}
\]
 represents a decaying oscillation along the $t$-axis, where the amplitude
decays exponentially. The function $\mathcal{L}[\phi_1(t)] = \Phi_1(s)$
has the asymptotic representation  $\Phi_1(s) \sim\frac{b}{c}s^{\alpha - 1}$
as $s\to 0$ and hence the function
$f_2(t) =  \int_0^\infty K_\alpha(\tau)e^{-\tau t}d\tau$ has the asymptotic
representation $f_2(t) \sim \frac{b}{c}\frac{t^{-\alpha}}{\Gamma(1-\alpha)}$
as $t \to \infty$. The derivative $u'(t)$ has a similar asymptotic
representation
$u'(t) \sim \frac{b}{c}\frac{(-\alpha) t^{-\alpha - 1}}{\Gamma(1-\alpha)}$ as
$t \to \infty$. We conclude that the solution $u(t)$ of the equation \eqref{FVE}
decays towards  0  as $t \to \infty$ like $t^{- \alpha }$  and $u'(t)$
has similar asymptotic properies. This means that the equilibrium
 $x = (x_1, x_2) = (u, u') = (0, 0)$ of the system of equations, corresponding
to the the solution $u(t)$ of the equation \eqref{FVE}, is asymptotically stable,
but not exponentially.

We were motivated by the paper \cite{SH}, where an existence and uniqueness
result for the initial value problem
 \begin{equation}\label{S}
 Au'' + \sum_{k=1}^NB_k\,^cD^{\alpha_k}u(t) = f(t),\quad
 u(0) = u_0,\quad u'(0) = c_1
 \end{equation}
with $0 < \alpha_k < 2$, $k = 1, 2,\dots, N$  is proved. The Caputo fractional
derivatives in the equation \eqref{S} play there the role of damping terms.

In the paper \cite{FT} the  initial value problem
\begin{equation}\label{IN}
^{RL}D^\alpha x(t) = f(t, x(t)),\quad t > 0,\quad
\lim_{t\to 0}t^{1-\alpha}x(t) = b,\quad \alpha \in (0, 1),\quad b \in \mathbb{R},
\end{equation}
where $f\colon \mathbb{R}^+ \times \mathbb{R}\to \mathbb{R}$,
$\phi\colon R^+ \to R^+$ are continuous functions, is studied.
It is assumed there that
\begin{equation}\label{PC}
    |f(t, x)| \leq t^\mu\phi(t)e^{-\sigma t}|u|^m,\quad
\text{for all }  (t, u) \in \mathbb{R}^+ \times \mathbb{R},
\end{equation}
where $\mu \geq 0$, $m > 1$, $\sigma > 0$. Using the desingularization method,
proposed in \cite{MED1} and the Bihari inequality, it is proven there that if
\[
\|\phi\|_q = \int_0^\infty\phi(s)^qds
< L:= \frac{\Gamma(\alpha)}{2^{1+m-\alpha }|b|^{m-1}}
\Big( \frac{2^m}{m-1} \Big)^{1/q}
\Big[\frac{(\sigma p)^{\lambda_1}}{\Gamma(\lambda_1)
(1 + \frac{\lambda_1}{\lambda_2})}\Big]^{1/p},
\]
where $pq = p + q$,  $\lambda_1 = 1+ p(\mu -(1 - \alpha)m]$,
$\lambda_2 = 1 + p(\alpha -1)$, then any solution $x(t)$ of the initial
value problem \eqref{IN} is global and there exists a constant $c >0$
such that $|x(t)| \leq \frac{c}{t^{1-\alpha}}$ for all $t \in (0, \infty)$.
This means that the trivial solution $x_0(t) \equiv 0$ is asymptotically stable.
  It is obvious that similar results for more general type of power
nonlinearities are extraordinary complicated.
     In the paper \cite{MED2} the equation
\begin{equation}\label{E2}
    \dot x(t) =Ax(t) + f\big(t, x(t), ^{RL}I^{\alpha_1}[g_1x](t),
\dots, ^{RL}I^{\alpha_m}[g_mx](t)\big),\quad x(t) \in  \mathbb{R}^N,
\end{equation}
where $ f\colon \mathbb{R}  \times \mathbb{R}^N  \times  \mathbb{R}^N
\to \mathbb{R}^N  $  is a continuous mapping,
$g_i\colon \mathbb{R}  \times \mathbb{R}^N \to \mathbb{R}^N ,
(t, x) \mapsto g_i(t, x)$,  $i=1,2,\dots,m$ are continuous mappings and
\begin{equation}
    ^{RL}I^{\alpha_i}[g_ix](t) := \frac{1}{\Gamma(\alpha_i)}
\int_0^t(t-s)^{\alpha_i-1}g_i(s, x(s))ds,\quad 0 < \alpha_i < 1,\; i=1,2,\dots, m
\end{equation}
is studied. A sufficient condition for the exponential stability of the
trivial solution $x(t) \equiv 0$ of this equation is proven there.
In the paper \cite{MED3} a sufficient condition for the non-existence of
blow-up solutions for a fractional functional-differential equations of the form
 \begin{equation}
\begin{aligned}
  &\dot x(t) = Ax(t) + h\Big(t, x(t), x_t, (I^{\alpha_1}[g_1x])(t),
    \dots,  (I^{\alpha_m}[g_mx])(t)\Big),\quad t > 0,\\
    &x(t) = \Phi(t),\quad t \in [- r, 0],
\end{aligned}
\end{equation}
where $ r > 0$, $\Phi \in C_r:= C([- r, 0], X)$, $X$ is a Banach space,
$x(t) \in X$, $x_t \in C$, $x_t(\Theta):= x(t + \Theta)$ $t > 0$,
$\Theta \in [- r, 0]$, $A$ is the
infinitesimal generator of a strongly continuous semigroup
$\{S(t)\}_{t \geq 0}$, $S(t) := e^{At}$,
$h\colon \mathbb{R}_+ \times X \times C_r \times X^m\to X$,
$X^m := X\times \dots \times X$
($m$\,times)  is a continuous map, $\mathbb{R}_+ = [0, \infty)$,
$g_i\colon \mathbb{R}_+ \times X \to X$, $(t, x) \mapsto g_i(t, x)$,
$i=1,2,\dots,m$ are continuous maps, is proved


In this paper, we study equation \eqref{E2} with $g_i(t, x) \equiv x(t)$, i.e.,
we have the Riemann-Liouville fractional integrals of $x(t)$ in the
equation \eqref{E4} instead of the nonlinear functions $g_i(t, x(t))$.
Moreover,  the mapping $f$ is more general than in \cite{MED2}.
The aim is to give some conditions under which the trivial solution of
this equation is exponentially stable.



\section{Existence and uniqueness result}\label{s:2}

In this section, we prove a local existence and uniqueness result concerning
the initial value problem
\begin{equation}\label{IVP}
\begin{gathered}
\dot x(t) =Ax(t) + f\Big(t, x(t), ^{RL}I^{\alpha_1}x(t), \dots,
 ^{RL}I^{\alpha_m}x(t)\Big), \quad t > 0,\\
x(t) \in  R^N,\quad x(t_0) = x_0.
\end{gathered}
\end{equation}
Many papers are devoted to the fractional initial value problem
\begin{equation}\label{IVP1}
^{RL}D^\alpha x(t) = f(t, x(t)),\quad x(t_0) = x_0,
\end{equation}
where the right-hand side is independent of fractional derivatives and
fractional integrals, respectivelly. Some basic existence results for the
problem \eqref{IVP1} can be found in the monograph \cite{LAK},
where they are proved by using the classical Picard method of successive
 approximations. Many local and global existence results for various
classes of fractional differential equations are proved by using fixed
point theorems (see, e.g., the monograph \cite{ZHO}, \cite{FT1} and the
paper \cite{ZHO1}). Some existence results for the second-order abstract
differential equations on Banach spaces, involving several fractional
derivatives on their right-hand sides are proved in the papers
\cite{KMT,TA1,TA2}. In its proof, we use the method of Picard successive
 approximations. There is a problem to apply the Banach fixed point
theorem without the assumption of the global boundedness of the mapping $f$,
 because there are fractional integrals of the unknown function $x(t)$ in
its arguments.

\begin{theorem}\label{T1}
Let $G \subset  \mathbb{R} \times \mathbb{R}^N$ be a region,
$H_m \subset \mathbb{R}^m$ is a region with $0 \in H_m$ and
$f \in C(G \times H_m, \mathbb{R}^N)$ be a continuous locally Lipschitz mapping.
Then for any $(t_0,x_0) \in G, t_0 \geq 0, $ there exists a $\delta > 0$
such that the initial value problem \eqref{IVP} has a unique solution $x(t)$
on the interval $I_\delta = [t_0, t_0+\delta)$.
\end{theorem}

\begin{proof}
Let
\begin{equation}
\begin{aligned}
    G_0 &= \big\{(t, x, u_1,\dots, u_m) \in G\times H_m:
     t_0 \leq t \leq t_0 + a, t_0 \geq 0,\\
&\quad  \|x-x_0\| \leq b,\,  \|u_i\| \leq \|x_0\| +b, \, i=1, 2,\dots, m\big\},
\end{aligned}
\end{equation}
for some $a > 0$, $b > 0$. Let
\[
M_1 = \max_{\|x-x_0\| \leq b}\|Ax\|,\quad
M_2 = \max_{(t, x, u_1,\dots, u_m) \in G_0}\|f(t, x, u_1,\dots, u_m)\|
\]
 and  the mapping $f$ satisfies the condition
\begin{equation}\label{LLC}
\|f(t, x, u_1, u_2,\dots, u_m) - f(t, y, v_1, v_2,\dots, v_m)\|
 \leq   L_0\|x-y\| + \sum_{i=1}^mL_i\|u_i-v_i\|
\end{equation}
for all $(t, x, u_1, u_2,\dots, u_m), (t, y, v_1, v_2,\dots, v_m) \in G_0$. Let
\[
 0 < \delta = \min \big\{a, \frac{b}{M_1 + M_2},\, c,\,
\frac{1}{\|A\|+ L_0 + \sum_{i=1}^mL_i}\big\},
\]
  where $ c = \min_{1\leq i \leq m}\big [ \Gamma(\alpha_i)\alpha_i
\big]^{\frac{1}{\alpha_i}}$. Let $C_\delta := C(I_\delta, \mathbb{R}^N)$
be the Banach space of continuous mappings from $I_\delta$ into
$\mathbb{R}^N$ endowed with the metric
$d(h, g) :=\|h-g\|:= \max_{t \in I_\delta}\|h(t) - g(t)\|$.
Let us define the successive approximations
 $\{x_n\}_{n=0}^\infty, x_n \in C_\delta := C(I_\delta, \mathbb{R}^N)$,
$I_\delta = [t_0, t_0+\delta]$, by
\begin{equation}\label{PA}
\begin{gathered}
 x_0(t) \equiv x_0,\\
\begin{aligned}
&x_{n+1}(t) \\
&=x_0 + \int_{t_0}^tAx_n(s)ds 
 +  \int_{t_0}^tf\Big(s, x_n(s), \frac{1}{\Gamma(\alpha_1)}\\
&\quad\times \int_0^s(s-\tau)^{\alpha_1-1}x_n(\tau)d\tau,
\dots,  \frac{1}{\Gamma(\alpha_m)}
 \int_0^s(s-\tau)^{\alpha_m-1}x_n(\tau)d\tau\Big)ds,\\
& n = 0, 1, 2,\dots \,\,t \in I_\delta.
\end{aligned}
\end{gathered}
\end{equation}

First, let us prove that $\|x_n(t) - x_0\| \leq b$ for all
$n \geq 1$, $t \in I_\delta$.
From the definition of the number $c$  it follows that
\begin{equation}\label{I}
 \frac{1}{\Gamma(\alpha_i)}\int_0^t(t-s)^{\alpha_i-1}ds
\leq \frac{1}{\Gamma(\alpha_i)}\frac{\delta^{\alpha_i}}{\alpha_i} \leq
\frac{1}{\Gamma(\alpha_i)}\frac{c^{\alpha_i}}{\delta^{\alpha_i}}
\leq  \frac{1}{\Gamma(\alpha_i)}\frac{\Gamma(\alpha_i)\alpha_i}{\alpha_i} = 1
\end{equation}
for $ i = 1, 2,\dots, m$;
and so, we have
\begin{equation}
 \|\frac{1}{\Gamma(\alpha_i)}\int_{t_0}^t(t-\tau)^{\alpha_i-1}x_0d\tau \|
\leq  \frac{1}{\Gamma(\alpha_i)}\frac{\delta^{\alpha_i}}{\alpha_i}[\|x_0\|+ b]
 \leq   \|x_0\|+ b,
\end{equation}
for $i = 1, 2,\dots, m$, $t \in I_\delta$, and $i=1,2,\dots m$.
Hence, the first approximation $x_1(t)$ is well defined and
\[
\|x_1(t) - x_0\| \leq M_1\delta + M_2\delta = (M_1 + M_2)\delta
\leq (M_1 + M_2)\frac{b}{M_1 + M_2} = b,
\]
for $t \in I_\delta$. This yields the inequality
\[
\|x_1(t)\| \leq\|x_0\| + b\quad \text{for all }t \in I_\delta.
\]
and thus
\[
\Big(t, x_1(t), \frac{1}{\Gamma(\alpha_1)}
\int_0^t(t-\tau)^{\alpha_1-1}x_1(\tau)d\tau,\dots,
 \frac{1}{\Gamma(\alpha_m)}\int_0^t(t-\tau)^{\alpha_m-1}x_1(\tau)d\tau\Big) \in G_0
\]
for all $t \in I_\delta$.
Now, we find by using the Lipschitz condition \eqref{LLC} and
inequality \eqref{I} that
\begin{equation}
\begin{aligned}
\|x_2(t) - x_1(t)\|
&\leq \delta (\|A\| + L_0)\|x_1(t) - x_0(t)\|\\
&\quad  +  \sum_{i=1}^m\frac{L_i}{\Gamma(\alpha_i)}
\int_{t_0}^t\int_0^s(s-\tau)^{\alpha_i-1}\|x_1(\tau)- x_0(\tau)\|d\tau ds\\
&\leq \delta (\|A\| + L_0)\|x_1 - x_0\| \\
&\quad  + \sum_{i=1}^m\frac{L_i}{\Gamma(\alpha_i)}
\Big (\int_{t_0}^t\int_0^s (s-\tau)^{\alpha_i -1}d\tau ds\Big )\|x_1-x_0\|\\
&\leq \delta k \|x_1-x_0\|,
\end{aligned}
\end{equation}
where $k = \|A\| + L_0 + \sum_{i=1}^mL_i$ and so, we get
\[
\|x_2-x_1\| \leq k\delta \|x_1-x_0\|.
\]
Now assume that the estimate
\[
\|x_n(t) - x_{n-1}(t)\| \leq (k\delta )^{n-1}
\]
holds for $n > 2$. Then, using this inequality, the Lipschitz condition
\eqref{LLC} and the inequality \eqref{I}, one can get
\[
\|x_{n+1}(t) - x_n(t)\| \leq (k\delta )^n\|x_1-x_0\|
\]
and so, we have
\[
\|x_{n+1} - x_n\| \leq (k\delta )^n\|x_1-x_0\|.
\]
Since
\begin{equation}
x_n(t) = x_0(t) + \sum_{i=1}^n[x_i(t) - x_{i-1}(t)]\ \ \text{with}\ \ x_0(t) \equiv x_0,
\end{equation}
we obtain
\begin{equation}
\begin{aligned}
\|x_0(t) + \sum_{i=1}^n[x_i(t) - x_{i-1}(t)]\|
&\leq   \|x_0\| + \sum_{i=1}^n\|x_i(t) - x_{i-1}(t)\|  \\
&\leq\Big ( \|x_0\| +  \sum_{i=1}^n(k\delta)^i\Big )\|x_1-x_0\|,
\end{aligned}
\end{equation}
for all $t \in I_\delta$.
From the definition of $\delta$ it follows that  $k\delta < 1$, and
so the series $\|x_0\| + \sum_{i=1}^\infty (k\delta)^i$ is convergent.
 This yields the uniform convergence of the sequence $\{x_n(t)\}_{i=0}^\infty$
on the interval $I_\delta$ to a continuous mapping $x \in C_\delta$.
This implies
\begin{equation}
\begin{aligned}
& \lim_{n \to \infty}\frac{1}{\Gamma(\alpha_i)}
 \int_0^s(s-\tau)^{\alpha_i-1}x_n(\tau)d\tau \\
& =\frac{1}{\Gamma(\alpha_i)}\int_0^s(s-\tau)^{\alpha_i-1}x(\tau)d\tau
\quad\text{for } i = 1, 2,\dots, m,\; s \in I_\delta
\end{aligned}
\end{equation}
and therefore
\begin{equation}
\begin{aligned}
&\lim_{n \to \infty} f\Big(s, x_n(s),
 \frac{1}{\Gamma(\alpha_1)}\int_0^s(s-\tau)^{\alpha_1-1}x_n(\tau)d\tau, \\
&\dots, \frac{1}{\Gamma(\alpha_m)}
 \int_0^s(s-\tau)^{\alpha_m-1}x_n(\tau)d\tau\Big)ds\\
&=  f\Big(s, x(s),\frac{1}{\Gamma(\alpha_1)}
 \int_0^s(s-\tau)^{\alpha_1-1}x(\tau)d\tau, \\
&\dots, \frac{1}{\Gamma(\alpha_m)}\int_0^s(s-\tau)^{\alpha_m-1}x(\tau)d\tau\Big)ds
\quad \text{for all } s \in I_\delta.
\end{aligned}
\end{equation}
 Therefore from \eqref{PA} it follows that $x(t)$ is a solution of the initial
value problem \eqref{IVP}, defined on the interval $I_\delta$.
Now let us prove its uniqueness. Assume that there are two different solutions
$x, y \in C_\delta$ of the initial value problem \eqref{IVP}.
 Let  $w(t) := \|x(t) - y(t)\|$, $t \in I_\delta$ and
$W = \max_{t \in I_\delta}w(t)$. Then, by using the Lipschitz condition \eqref{LLC}
and the inequality \eqref{I} we obtain
\begin{equation}
\begin{aligned}
    w(t)
&\leq (\|A\|+ L_0 )\int_{t_0}^tw(s)ds
 + \sum_{i=0}^m\frac{L_i}{\Gamma(\alpha_i)}
 \int_{t_0}^t\int_0^s (s-\tau)^{\alpha_i-1}\|x(\tau) - y(\tau)\|d\tau ds\\
&\quad +\delta \Big(\|A\|+ L_0+  \sum_{i=0}^m\frac{L_i}{\Gamma(\alpha_i)}
\int_{t_0}^t(t-\tau)^{\alpha_i-1}d\tau \Big)W\\
&\leq \delta \Big (\|A\|+ L_0+ \sum_{i=0}^m L_i\Big )W \\
&= (\delta k)W\quad\text{for all } t \in I_\delta
\end{aligned}
\end{equation}
and this yields the inequality $W \leq (k \delta) W < W$.
This is a contradiction and hence, we have $x(t) =y(t)$ for all $t \in I_\delta$.
\end{proof}


\section{Stability Theorem}\label{s:3}

In this section, we prove a result on the exponential stability of the
trivial solution $x(t) \equiv 0$ of the equation \eqref{E4}. In its proof,
we apply a desingularization method, proposed in the paper \cite{MED1},
where it is applied in the study of nonlinear integral inequalities
with weakly singular kernels.

We assume that the following conditions are satisfied:
\begin{itemize}
\item[(A1)]
 \[
\|e^{At}x\| \leq Ke^{-at}\|x\|\quad \text{for all }
 t \geq 0,\ x \in \mathbb{R}^N,
\]
where $K > 0$, $a > 0$ are constants;

\item[(A2)] The mapping $f\colon \mathbb{R}\times
\mathbb{R}^N \times \mathbb{R}^{mN}  \to \mathbb{R}^N $ is continuous and
it satisfies  the condition
\begin{equation}
\begin{aligned}
    \|f(t,  x, u_1, u_2,\dots, u_m)\|
&\leq  \mu_1(t)e^{-\gamma_1 t}\|x\|+ \sum_{i=1}^m\lambda_{i1}(t)
 e^{-\gamma_{i1}t}\|u_i\| \\
&\quad +     \sum_{j=2}^l\mu_j(t)\|x\|^{k_j} + \sum_{j=2}^l
\sum_{i=1}^m\lambda_{ij}(t)\|u_i\|^{k_j}
\end{aligned}
\end{equation}
for all $(t, x, u_1, u_2,\dots, u_m) \in  \mathbb{R} \times \mathbb{R}^N
\times \mathbb{R}^{mN}$, where $\mu_j(t), \lambda_{ij}(t)$,
$i=1, 2,\dots, m$, $j=1,2,\dots, l$ are nonnegative continuous functions on
$[0,\infty)$, $\gamma_1 > 0$,  $\gamma_{i1} > a+1$,
$i = 1, 2,\dots, m$,  $1 = k_1 < k_2 <\dots <k_l$,
$\|z\| = \max\{|z_1|,|z_2|,\dots, |z_N|\}$;

  \item[(A3)]
 There exist numbers $p_i > 1$,  $i=1,2,\dots,m$ such that
 \[
    p_i(\alpha_i-1) + 1 >0, \quad i =1, 2,\dots, m
\]
and
\begin{equation}
\omega_j := \int_0^\infty \mu_j(s)^qds < \infty,\quad j =1, 2, 3, \dots, l,
\end{equation}
where $q=q_1q_2\dots q_m$,   $q_i = \frac{p_i}{p_i-1}$,  $i=1, 2, \dots, m$;

\item[(A4)]
\begin{equation}
\begin{gathered}
\eta_{i1} := \int_0^\infty e^{-[\gamma_{i1} - (a+1)]s}\lambda_{i1}(s)ds < \infty, \\
\eta_{ij} :=\int_0^\infty e^{(a+k_j)s}s^{\frac{k_j-1}{q_i}}\lambda_{ij}(s)ds
< \infty,\quad i=1, 2, \dots, m,\; j =  2, 3, \dots, l,
\end{gathered}
\end{equation}
where $q_1, q_2,\dots, q_m$ are defined as in (A3).

\item[(A5)]
The mapping $f(t,  x, u_1, u_2,\dots, u_m)$ is locally lipschitz with
respect to the variables $x, u_1,\dots, u_m$.
\end{itemize}

In the proof of the main result, we use the following corollary of the
 Pinto's inequality (see  \cite[Theorem 1]{PIN}, \cite[Theorem 10.2]{BA}
and \cite[Example 5]{TA}).  We present it in the form of the next lemma
also with its proof, because we did not find this formulation in literature.

\begin{lemma}\label{L1}
Let $c >0$ be a constant, $\Psi_j(t)$, $j = 1, 2,\dots, l$ be continuous,
 nonnegative functions on $[a, \infty)$ and $u(t)$ be a continuous
nonnegative function satisfying the integral inequality
\[
u(t) \leq c + \sum_{j=1}^l\int_a^t\Psi_j(s)u(s)^{k_j}ds,\quad t \in [a, \infty),
\]
where $a \in \mathbb{R}$, $1 = k_1 \leq k_2 <\dots \leq k_l$.
Let the following conditions be satisfied:
\begin{equation}\label{D}
 (k_j-1)(cD_j)^{k_j-1}\int_a^\infty\Psi_js)ds < 1,\quad j=2, 3, \dots, l,\quad
 \int_0^\infty \Psi_1(s)ds < \infty,
\end{equation}
where
\begin{equation}\label{CAS}
\begin{gathered}
    D_1 =e^{\int_0^\infty\Psi_1(s)ds},\\
    D_j= \Big(1-(k_j-1)(D_{j-1}c)^{k_j-1}\int_a^\infty\Psi_j(s)ds
\Big)^{-\frac{1}{k_j-1}},\quad j = 2, 3,\dots, l.
\end{gathered}
\end{equation}
Then
\begin{equation}\label{PI}
u(t) \leq cD_1D_2\dots D_l,\quad\text{for all } t \in [a,\infty).
\end{equation}
\end{lemma}

\begin{proof}
By \cite[Theorem 1]{PIN} (see also \cite[Example 5]{TA}),
\begin{equation}\label{PIN}
u(t) \leq W_l^{-1}\Big[W_l(c_{l-1}(t))+ \int_a^t\Psi_l(s)ds\Big],
\end{equation}
where
\begin{equation}
\begin{gathered}
 c_0(t)\equiv c,\quad c_i(t) = W_i^{-1}\Big[W_i(c_{i-1}(t)) 
+ \int_a^t\Psi_i(s)ds\Big], \\
 W_i(z) = \int_{u_i}^z\frac{dy}{y^{k_i}},\quad  y \geq u_i > 0,\quad
 i = 1, 2,\dots ,l.
\end{gathered}
\end{equation}
One can calculate that
\[
W_1(y) = \frac{1}{1-k_1}\big[y^{1-k_1} - c^{1-k_1}\big],\quad  
W_1^{-1}(u) = \big[ c^{1-k_1} - (k_1-1)u\big]^{-\frac{1}{k_1-1}}
\]
and
\begin{equation}\label{D1}
c_1(t) = W_1^{-1}\Big[W_1(c) + \int_a^t\Psi_1(s)ds\Big] \leq cD_1.
\end{equation}
From the assumption \eqref{D} it follows that $0 < D_1 < \infty$. 
Using the inequality \eqref{D1}, we obtain
\begin{equation}\label{D2}
\begin{aligned}
c_2(t) 
&\leq W_2^{-1}\big [W_2(c_1(t))+ \int_a^t\Psi_2(s)ds\big ] 
 \leq W_2^{-1}\big [W_2(cD_1) + \int_a^t\Psi_2(s)ds\big ]\\
& \leq \Big[ (D_1c)^{1-k_2} - (k_2-1)\int_a^t\Psi_2(s)ds 
\Big]^{-\frac{1}{k_2-1}} \leq cD_1D_2,
\end{aligned}
\end{equation}
where
\[
D_2 = \Big[ 1 - (k_2-1)(cD_1)^{k_2-1}\int_a^\infty\Psi_2(s)ds 
\Big]^{-\frac{1}{k_2-1}}.
\]
Now, assume that
\[
c _{l-1}(t) \leq cD_1D_2\cdots D_{l-1}.
\]
Using the same arguments as above one can prove inequality \eqref{PI}.
\end{proof}

\begin{theorem}\label{T}
Let  conditions {\rm (A1)--(A5)} be satisfied and $\|x_0\| < \rho$, where 
$\rho = \infty$, if $l = 1$  and if $l>1$, then
\begin{equation}\label{RHO}
    \rho = \sup\big\{z \in \mathbb{R}: |C(z)D_{j-1}(z)| 
< \Big[\frac{1}{(k_j-1)G_j}\Big ]^{\frac{1}{k_j-1}},\; i = 2, 3,\dots, l\big\},
\end{equation}
where
\begin{equation}
\begin{gathered}
 C(z) = d^qK^qz^q,\quad d = m(l+1) + 2,\quad D_1(z) \equiv G_1,\\
 D_2(z) = \big[1 - (k_2 - 1)[C(z)D_1(z)]^{k_2-1}G_2\big]^{- \frac{1}{k_2-1}},\\
 D_j(z) = \big[1 - (k_j - 1)\big [C(z)D_{j-1}(z)\big ]^{k_j-1} G_j
 \big]^{- \frac{1}{k_j-1}}, \quad  j = 3,\dots, l,
\end{gathered}
\end{equation}
\begin{equation}
\begin{gathered}
  G_1= d^{q-1}K^q\Big\{\Big(\frac{1}{\gamma_1p}\Big)^{q/p}\omega_1 +
    \sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^q\eta_{i1}^q
 \Big(\frac{1}{aq_i\hat{p}_i}\Big)^{\frac{1}{\hat{p}_i}}\frac{1}{q}\Big\},
\\
  G_j= d^{q-1}K^q\Big\{\Big[\frac{1}{(k_j-1)ap}\Big]^{q/p}\omega_j +
    \sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^{k_jq}\eta_{ij}^q
\Big (\frac{1}{aq_i\hat{p}_ik_j}\Big)^{\frac{1}{\hat{p}_i}}\frac{1}{k_jq}\Big\},
\end{gathered}
\end{equation}
$j = 2, 3, \dots , l$, $p_i > 1$, $q_i = \frac{p_i}{p_i - 1}$, $i=1,2,\dots, m$,
$q = q_1q_2\dots q_m$, $p = \frac{q}{q-1}$, 
${\hat q}_i = q_1\dots q_{i-1}q_{i+1}\dots q_m$, 
$\hat{p}_i = \frac{{\hat q}_i }{{\hat q}_i -1}$,
\begin{equation}\label{Q}
Q_i = \Big(\frac{\Gamma(p_i(\alpha_i-1)+1)}{p_i^{p_i(\alpha_i-1)+1}}
\Big)^{1/p_i}.
\end{equation}
Then for the solution $x(t)$ of the equation \eqref{E4}, satisfying the 
condition $x(0) = x_0$, the following inequality holds:
\begin{equation}
\|x(t)\| \leq \Psi(\|x_0\|)e^{-at}\quad\text{for all } t \in [0, \infty),
\end{equation}
where the function $\Psi(w)$ is defined as
\[
\Psi(w) = \Big[C(w)D_1(w)D_2(w)D_2(w)\dots D_l(w)\Big]^{1/q},\quad |w| < \rho,
\]
for which $\lim_{w \to 0}\Psi(w) = 0$ and if $l > 1$, then 
$\lim_{w \to\rho^-}\Psi(w) = \infty$.
\end{theorem}

\begin{proof} 
Let $x(t)$ be the maximal solution of the equation \eqref{E4}, 
defined on the interval $[0, d)$, 
satisfying the condition $x(0) = x_0 \in \mathbb{R}^N$. 
From Theorem \ref{T1} it follows that this solution exists. Then
\[
x(t) = e^{At}x_0 + \int_0^te^{A(t-s)}f\big (s, x(s), ^{RL}I^{\alpha_1}x(s),
\dots, ^{RL}I^{\alpha_m}x(s)\big )ds,\quad t \in [0, d)
\]
and the conditions (A1), (A2) yield
\begin{equation}
\begin{aligned}
    \|x(t)\|
&\leq Ke^{-at}\|x_0\|  + Ke^{-at}\int_0^te^{-(\gamma_1 - a)s}\mu_1(s)\|x(s)\|ds  \\
&\quad + Ke^{-at}\sum_{i=1}^m\int_0^te^{-(\gamma_{i1} a)s}
 \lambda_{i1}(s)\|^{RL}I^{\alpha_i}x(s)\|ds    \\
&\quad +Ke^{-at}\int_0^t e^{as}\Big(\sum_{j=2}^l\mu_j(s)\|x(s)\|^{k_j}
 + \sum_{j=2}^l\sum_{i=1}^m\lambda_{ij}(s)\|^{RL}I^{\alpha_i}x(s)\|^{k_j} \Big)ds.
\end{aligned}
\end{equation}
If $u(t) = e^{at}\|x(t)\|$, then  $\|x(t)\| = e^{-at}u(t)$,
$\|x(t)\|^{k_j} = e^{-ak_jt}u(t)^{k_j}$,
\begin{equation}
\begin{aligned}
    \|^{RL}I^{\alpha_i}x(s)\|
&=\frac{1}{\Gamma(\alpha_i)}  \Big\|\int_0^t(t-s)^{\alpha_i-1}x(s)ds\Big \| \\
&\leq   \frac{1}{\Gamma(\alpha_i)}\int_0^t(t-s)^{\alpha_i-1}\|x(s)\|ds\\
&\leq \frac{1}{\Gamma(\alpha_i)}\int_0^t(t-s)^{\alpha_i-1}e^{-as}u(s)ds,
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
 \|^{RL}I^{\alpha_i}x(s)\|^{k_j}
&= \frac{1}{\Gamma(\alpha_i)^{k_j}}\Big \|\int_0^t(t-s)^{\alpha_i-1}
 x(s)ds\Big\|^{k_j} \\
& \leq  \frac{1}{\Gamma(\alpha_i)^{k_j}}\Big( \int_0^t(t-s)^{\alpha_i-1}\|x(s)\|ds
 \Big)^{k_j}\\
&\leq \frac{1}{\Gamma(\alpha_i)^{k_j}}\Big( \int_0^t(t-s)^{\alpha_i-1}
 e^{-as}u(s)ds \Big)^{k_j}.
\end{aligned}
\end{equation}
Now, we apply the desingularization method as follows:
\begin{equation}
\begin{aligned}
& \int_0^t(t-s)^{\alpha_i-1}e^se^{-s}e^{-as}u(s)ds  \\
&\leq\Big( \int_0^t(t-s)^{p_i(\alpha_i-1)}e^{p_is}ds  \Big)^{1/p_i}
\Big (\int_0^te^{-q_is}e^{-aq_is}u(s)^{q_i}ds\Big)^{1/q_i}
\end{aligned}
\end{equation}
with $q_i = \frac{p_i}{p_i-1}$, $p_i(\alpha_i-1)+1 > 0$.
By the following inequality, proved in \cite{MED1} (see also \cite{MED4}),
\[
\Big(\int_0^t(t-s)^{p_i(\alpha_i-1)}e^{p_is}ds\Big)^{1/p_i} \leq Q_ie^t,
\]
where the number $Q_i$ is given by \eqref{Q},we obtain the inequality
\begin{equation}
\Big( \int_0^t(t-s)^{\alpha_i-1}\|x(s)\|ds \Big)^{k_j}
\leq Q_i^{k_j}e^{k_jt}\Big (\int_0^te^{-q_is}e^{-aq_is}u(s)^{q_i}ds
\Big)^{k_j/q_i}.
\end{equation}
Then the H\"older inequality yields
\[
\Big(\int_0^te^{-q_is}e^{-aq_is}u(s)^{q_i}ds\Big)^{k_j}
\leq t^{k_j-1}\int_0^te^{-q_ik_js}e^{-aq_ik_js}u(s)^{k_jq_i}ds,\quad j > 2
\]
and hence, we have the inequality
\[
 \|^{RL}I^{\alpha_i}x(s)\|^{k_j}
\leq  \Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^{k_j}e^{k_js}
s^{\frac{k_j-1}{q_i}}\Big(\int_0^se^{-q_ik_j\tau}e^{-aq_ik_j\tau}
u(\tau)^{k_jq_i}d\tau\Big )^{1/q_i}.
\]

Thus, we obtain the inequality
\begin{align*}
u(t) &\leq K\|x_0\| + K\int_0^t e^{-\gamma_1s}\mu_1(s)u(s)ds \\
    &\quad + K \sum_{i=1}^m\frac{Q_i}{\Gamma(\alpha_i)}
\int_0^te^{-[\gamma_{i1} - (a+1))s]}\lambda_{i1}(s)
\Big(\int_0^se^{-q_i\tau} e^{-aq_i\tau}u(\tau)^{q_i}d\tau\Big )^{1/q_i}ds\\
&\quad + K\sum_{j=2}^l\int_0^te^{-(k_j -1)a)s}\mu_j(s)u(s)^{k_j}\Big)ds
 + K\sum_{j=2}^l\sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^{k_j} \\
&\quad\times    \int_0^t\lambda_{ij}(s)e^{(a+k_j)s}s^{\frac{k_j-1}{q_i}}
\Big (\int_0^se^{-q_ik_j\tau}e^{-aq_ik_j\tau}u(\tau)^{k_jq_i}d\tau\Big)^{1/q_i}ds.
\end{align*}

If $q = q_1q_2\dots q_m$ and $d = m(l+1) + 2$, then the inequality 
$(z_1+z_2+\dots +z_d)^q \leq d^{q-1}(z_1^q+z_2^q+\dots +z_m^q)$,
 valid for any $z_1, z_2,\dots, z_d \geq 0$, yields
\begin{align*}
&u(t)^q  \\
&\leq d^{q-1}K^q\|x_0\|^q + d^{q-1}K^q\Big(\int_0^te^{-\gamma_1s}\mu_1(s)u(s)ds
  \Big)^q
\\
&\quad+ d^{q-1}K^q\sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^q
 \Big[\int_0^te^{-[\gamma_{i1}- (a+1)] s}\lambda_{i1}(s)
 \Big(\int_0^se^{-q_i\tau}e^{-aq_i\tau}u(\tau^{q_i}d\tau\Big)ds  \Big ]^q
\\
&\quad +d^{q-1}K^q\sum_{j=2}^l\Big[\int_0^te^{-(k_j-1)as}\mu_j(s)u(s)^{k_j}ds\Big]^q 
 + d^{q-1}K^q\sum_{j=2}^l\sum_{i=1}^m \Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^{k_jq}
\\
&\quad\times\Big[\int_0^t\lambda_{ij}(s)e^{(a+k_j)s}s^{\frac{k_j-1}{q_i}}
\Big (\int_0^se^{-aq_ik_j\tau}e^{-k_jq_i\tau}u(\tau)^{k_jq_i}d\tau\Big)^{1/q_i} 
ds\Big]^q \\
&\leq d^{q-1}K^q\|x_0\|^q +  d^{q-1}K^q\Big(\frac{1}{\gamma_1p}\Big)^{q/p}
\int_0^t\mu_1(s)^qu(s)^qds\\
&\quad + d^{q-1}K^q\sum_{i=1}^m\Big( \frac{Q_i}{\Gamma(\alpha_i)}\Big)^q
 \eta_{i1}\Big (\int_0^te^{-q_i\tau} e^{-aq_i\tau}u(\tau)^{q_i}d\tau 
 \Big)^{\hat{q}_i}\\
&\quad +d^{q-1}K^q\sum_{j=2}^l\sum_{i=1}^m \Big(\frac{Q_i}{\Gamma(\alpha_i)}
\Big)^{k_jq}\eta_{ij}^q
\Big(\int_0^te^{-q_ik_j\tau} e^{-aq_ik_j\tau}u(\tau)^{q_ik_j}d\tau \Big)^{\hat{q}_i},
\end{align*}
where $\hat{q}_i = q_1q_2\dots q_{i-1}q_{i+1}\dots q_m$.

Using  H\"older's inequality with 
$\hat{q}_i , \hat{p}_i = \frac{\hat{q}_i}{\hat{q}_i-1}$ and with $p, q$, 
 we obtain
\begin{equation}
\begin{aligned}
    \Big(\int_0^te^{-\gamma_1s}\mu_1(s)u(s)ds \Big)^q
    &\leq \Big(\int_0^te^{-\gamma_1ps}ds\Big )^{q/p}\int_0^t\mu_1(s)^qu(s)^qds
\\
    &\leq  \Big(\frac{1}{\gamma_1p}\Big)^{q/p}\int_0^t\mu_1(s)^qu(s)^qds,
\end{aligned}
\end{equation}
\[
 \Big(\int_0^te^{-[\gamma_{i1} - (a+1)]s}\lambda_{i1}(s)u(s)ds  \Big)^q
    \leq \Big (\frac{1}{[\gamma_{i1} - (a+1)]p}\Big)^{q/p}
 \int_0^t\mu_{i1}(s)^qu(s)^qds,
\]
\begin{align*}
     \Big(\int_0^te^{as}\mu_j(s)e^{-k_jas}u(s)^{k_j}ds\Big)^q 
&\leq \Big(e^{-(k_j-1)aps}ds\Big )^{q/p}\int_0^t\mu_j(s)^qu(s)^{k_jq}ds\\
& \leq\Big[\frac{1}{(k_j-1)ap}\Big]^{q/p}\int_0^t\mu_j(s)^qu(s)^{k_jq}ds,
\end{align*}
\begin{equation}
\begin{aligned}
& \Big[\int_0^t\lambda_{ij}(s)e^{(a+k_j)s}s^{\frac{k_j-1}{q_i}}
\Big (\int_0^se^{-aq_ik_j\tau}e^{-k_jq_i\tau}u(\tau)^{k_jq_i}d\tau\Big)^{1/q_i}
 ds\Big]^q\\
&\leq \Big(\int_0^t \lambda_{ij}(s)e^{(a+k_j)s}s^{\frac{k_j-1}{q_i}}ds \Big)^q
    \Big(\int_0^t e^{-aq_ik_j\tau}e^{-k_jq_i\tau}u(\tau)^{k_jq_i}d\tau 
\Big)^{\hat{q}_i},
 \end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
&\Big(\int_0^t e^{-aq_ik\tau}e^{-k_jq_i\tau}u(\tau)^{kq_i}d\tau \Big)^{\hat{q}_i}\\
&\leq  \Big(\int_0^te^{-aq_i\hat{p}_ik\tau}d\tau  \Big)^{\frac{1}{\hat{p}_i}}
\Big (\int_0^te^{-k_jq\tau}u(\tau)^{k_jq}d\tau\Big)\\
& \leq \Big(\frac{1}{aq_i\hat{p}_ik_j}\Big )^{\frac{1}{\hat{p}_i}}\int_0^t
    e^{-k_jq\tau}u(\tau)^{k_jq}d\tau.
 \end{aligned}
\end{equation}
The above inequalities yield
\begin{equation}
\begin{aligned}
&u(t)^q \\
&\leq d^{q-1}K^q \|x_0\|^q + d^{q-1}K^q\Big(\frac{1}{\gamma_1p}\Big )
\frac{q}{p}\int_0^t\mu_1(s)^qu(s)^qds \\
&\quad + d^{q-1}K^q\sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}
 \Big)^q\eta_{i1}^q\Big (\frac{1}{aq_i\hat{p}_i}\Big)^{\frac{1}{\hat{p}_i}}
 \int_0^te^{-q\tau}u(s)^qds\\
&\quad +  d^{q-1}K^q\sum_{j=2}^l\Big[\frac{1}{(k_j-1)ap}\Big ]^{q/p}
 \int_0^t\mu_j(s)^qu(s)^{k_jq}ds\\
&\quad  + d^{q-1}K^q\sum_{j=2}^l\sum_{i=1}^m
 \Big (\frac{Q_i}{\Gamma(\alpha_i)}\Big)^{k_jq}\eta_{ij}^q
 \Big (\frac{1}{aq_i\hat{p}_ik_j}\Big)^{\frac{1}{\hat{p}_i}}
 \int_0^te^{-k_jq\tau}u(\tau)^{k_jq}d\tau.
\end{aligned}
\end{equation}
Therefore   $v(t) = u(t)^q$ satisfies the integral inequality
\begin{equation}\label{G}
v(t) \leq d^{q-1}K^q\|x_0\|^q +\sum_{j=1}^l \int_0^tF_j(s)v(s)^{k_j}ds,
\end{equation}
where
\begin{gather*}
F_1(t) = d^{q-1}K^q\Big\{\Big(\frac{1}{\gamma_1p}\Big)^{q/p}\mu_1(t)^q 
+ \sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^q\eta_{i1}^q
\Big(\frac{1}{aq_i\hat{p}_i}\Big )^{\frac{1}{\hat{p}_i}}e^{-qt}\Big\},
\\
F_j(t) = d^{q-1}K^q\Big\{\Big[\frac{1}{(k_j-1)ap}\Big]^{q/p}\mu_j(s)^q 
 +    \sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^{k_jq}\eta_{ij}^q
\Big (\frac{1}{aq_i\hat{p}_ik_j}\Big)^{\frac{1}{\hat{p}_i}}e^{-k_jq\tau}\Big\}.
\end{gather*}
Obviously,
\begin{gather*}
\int_0^\infty F_1(s)ds=d^{q-1}K^q\Big\{\Big(\frac{1}{\gamma_1p}\Big)^{q/p}\omega_1
 +    \sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^q\eta_{i1}^q
\Big(\frac{1}{aq_i\hat{p}_i}\Big )^{\frac{1}{\hat{p}_i}}\frac{1}{q}\Big\},
\\
\int_0^\infty F_j(s)ds=d^{q-1}K^q\Big\{\Big[\frac{1}{(k_j-1)ap}\Big]^{q/p}\omega_j 
 + \sum_{i=1}^m\Big(\frac{Q_i}{\Gamma(\alpha_i)}\Big)^{k_jq}\eta_{ij}^q
\Big (\frac{1}{aq_i\hat{p}_ik_j}\Big)^{\frac{1}{\hat{p}_i}}\frac{1}{k_jq}\Big\},
\end{gather*}
i.e., $G_j= \int_0^\infty F_j(s)ds  < \infty$, $j = 1, 2,\dots, l$.

From Lemma \ref{L1} it follows that if $\|x_0\| < \rho$, where $\rho > 0$ 
is defined by \eqref{RHO}, then
\begin{equation}\label{V}
v(t) \leq \Psi_0\|x_0\|) 
= C(\|x_0\|)D_1(\|x_0\|)D_2(\|x_0\|)D_2(\|x_0\|)\dots D_l(\|x_0\|),
\end{equation}
for $t \in [0, d)$, where
\begin{equation}\label{UE}
\begin{gathered}
    C(\|x_0\|) = d^qK^q\|x_0\|^q, D_1(\|x_0\|) = \int_0^\infty F_1(s)ds,
\\
    D_2(\|x_0\|) = \Big[1 - (k_1 - 1)[C(\|x_0\|D_1(\|x_0\|)]^{k_1-1}
\int_0^\infty F_2(s)ds \Big ]^{- \frac{1}{k_1-1}},
\\
    D_j(\|x_0\|) = \Big[1 - (k_j - 1)\big [C(\|x_0\|)D_{j-1}(\|x_0\|)\big ]^{k_j-1} 
\int_0^\infty F_j(s)ds\Big]^{- \frac{1}{k_j-1}},
\end{gathered}
\end{equation}
for $ j = 3,\dots, l$.

Obviously, the right hand side of the inequality \eqref{V} is finite if 
$\|x_0\| < \rho$ and it is going to $\infty$ for $\|x_0\| \to \rho$, if $l > 1$. If $l= 1$, then it is defined for all $x_0 \in \mathbb{R}^N$.
 Hence, we have the estimate
\[
v(t) = u(t)^q =\big ( \|x(t)\|e^{at} \big )^q \leq \Psi_0(\|x_0\|),\, t \in [0, d),
\]
i.e.,
\begin{equation}\label{P}
\|x(t)\| \leq \Psi(\|x_0\|)e^{-at},\, t \in [0, d).
\end{equation}
From this inequality it follows that $\lim_{t \to d}x(t) = d_- < \infty$ and by Theorem \ref{T1} there is an $\epsilon > 0$ such that the inial value problem \eqref{IVP}
has a unique solution $w(t)$, defined on the interval $[0, d+\epsilon)$. This is a contradiction with the maximality of the solution $x(t)$. Hence, the inequality  \eqref{P} holds for all $t \in [0, \infty)$ and the proof is complete.
\end{proof}

\section{Illustrative example}\label{s:4}

Let us apply Theorem \ref{T} to the fractionally perturbed pendulum equation
\begin{equation}\label{EX}
    v''(t) + 2v'(t) + 4v(t) + 3e^{-4t}{^C}D^{1/3}v(t)
+ 5e^{-4t}\Big[ {^C}D^{1/2}v(t)\Big]^2  = 0.
\end{equation}
We can write this equation in the form of system \eqref{E4},
 where $x(t) = (x_1(t), x_2(t)) = (v(t), v'(t))$, $m =2$, 
$\alpha_1=1-\frac{1}{3} = \frac{2}{3}$, $\alpha_2 = 1 - \frac{1}{2} = \frac{1}{2}$,
\[
    A = \begin{pmatrix}
    0 & 1\\
    -4 & - 2
    \end{pmatrix},
\quad
    f(t, x, u_1, u_2) =\begin{pmatrix}
    0\\ -3e^{-4t}u_{12}  - 5e^{-4t}u_{22}^2,
    \end{pmatrix},
\]
where $u_1 = (u_{11},u_{12})$,  $u_2 =  (u_{21},u_{22})$,
\[
  f(t, x(t), {}^{RL}I^{\alpha_1}x(t), {}^{RL}I^{\alpha_2}x(t)) 
=\begin{pmatrix}
    0\\- 3e^{-4t}{^{RL}}I^{\frac{2}{3}}x_2(t) 
- 5e^{-4t}\Big[{^{RL}}I^{1/2}x_2(t)\Big]^2.
\end{pmatrix},
\]
Obviously,
\begin{equation}\label{If}
\|f(t, x, u_1, u_2)\| \leq 3e^{-4t}\|u_1\| +  5e^{-4t}\|u_2\|^2.
\end{equation}
The system has the form
\begin{equation}\label{ES}
    \dot x(t) = Ax(t) +
\begin{pmatrix}
    0\\
- 3e^{-4t}{^{RL}}I^{\frac{2}{3}}x_2(t) 
- 5e^{-4t}\Big[{^{RL}}I^{1/2}x_2(t)\Big]^2
\end{pmatrix}.
\end{equation}
Now, let us find the constants $a > 0$ and $K > 0$ from 
 condition (A1). The Jordan block of the matrix $A$ has the form
\[
\begin{pmatrix}
    \alpha & \beta\\ 
- \beta & \alpha,
    \end{pmatrix},
\]
where $\alpha =  - 1$, $\beta = 1$. One can prove by using the Putzer's 
method (see \cite{PUT}) that
\begin{equation}\label{ES1}
    e^{At}x = e^{\alpha t}\Big(\cos(\beta t) I 
+ \frac{1}{\beta}\sin (\beta t)(A- I)\Big)x,
\end{equation}
where $I$ is the unit matrix. This yields the estimate
\begin{equation}
    \|e^{At}x\| \leq e^{-t}\big (\|I\| + \|A-I\|\big )\|x\|,
\end{equation}
where we use the norm $\|C\| = \max \{|c_{11}| + |c_{12}|, |c_{21}| + |c_{22}|\}$ 
of a $2 \times 2$ matrix $C = (c_{ij})$. For this norm the inequality 
$\|Cy\| \leq \|C\|\|y\|$ is valid for any $y = (y_1, y_2)$ with the norm 
$\|y\| = \max \{|y_1|, |y_2|\}$. Since $\|I\| = 1, \|A-I\| = 7$, 
from \eqref{ES1} we obtain 
\begin{equation}
\|e^{At}x\| \leq 8 e^{-t}\|x\| \quad \text{for all } x \in \mathbb{R}^2,
\end{equation}
i.e., we have $a = 1$ and $K = 8$. This means that the condition (A1)
 is satisfied with $a = 1$,$ K = 8$.

From \eqref{If} it follows that the mapping $f$ has the form
\begin{equation}
\|f(t, x, u_1, u_2)\| \leq \lambda_{11}(t)e^{-\gamma_{11}t}\|u_1\|
 + \lambda_{12}(t)\|u_2\|^{k_2}
\end{equation}
with $\lambda_{11}(t) = 3$, $\gamma_{11} = 4$, $\lambda_{12}(t) = e^{-4t}, k_2 = 2$. 
This means that the condition (A2) is satisfied . The condition (A3)
is trivially satisfied because all $\omega_j$ are equal zero.
 The condition (A4) is also satisfied because
\begin{gather}\label{ETA1}
    \eta_{11}= \int_0^\infty  e^{-[\gamma_{11} 
- (a+1)]s}\lambda_{11}(s)ds = 3\int_0^\infty  e^{-2s}ds = \frac{3}{2},
\\
\label{ETA2}
\begin{aligned}
 \eta_{12}&= \int_0^\infty e^{[a+k_2]s}s^{\frac{k_2-1}{q}}\lambda_{12}(s)ds 
 = \int_0^\infty  e^{3s}s^{1/q}e^{-4s}ds \\
&= \int_0^\infty  s^{1/q}e^{-s} = \Gamma\Big(1 + \frac{1}{q}\Big)
\end{aligned}
\end{gather}
with $q = q_1q_2$, where we define $q_1 = 2$, $q_2 = 3$, i.e., $q = 6$. 
For the numbers $p_1 = \frac{q_1}{q_1 - 1}= 2$ and 
$p_2 = \frac{q_2}{q_2 - 1} = \frac{3}{2}$, $\hat{q}_1 = q_2 = 3$,
 $\hat{q}_2 = q_1 = 2, \hat{p}_1 = \frac{\hat{q}_1}{\hat{q}_1 - 1} 
= \frac{3}{2}$, $\hat{p}_2 = \frac{\hat{q}_2}{\hat{q}_2 - 1} = 2$,  we have
\[
p_1(\alpha_1 -1)+1 = \frac{1}{3},\quad p_2(\alpha_2 -1)+1 = \frac{1}{4}.
\]
Since the mapping $f$ is smooth in all its variables, the condition (A5)
 follows from the Lagrange mean value theorem.

Now let us calculate the numbers $Q_1, Q_2, G_1, G_2$ from the assumptions 
of Theorem \ref{T}:
\begin{gather}\label{Q1}
    Q_1 = \Big[\frac{\Gamma(p_1(\alpha_1 -1)+1)}{p_1^{p_1(\alpha_1 -1)+1}} 
\Big]^{1/p_1} 
= \Big[\frac{\Gamma(\frac{1}{3})}{2^{1/3}}  \Big]^{1/2}
 = \frac{\Gamma(\frac{1}{3})^{1/2}}{2^{1/6}}, \\
\label{Q2}
    Q_2 = \Big[\frac{\Gamma(p_2(\alpha_2 -1)+1)}{p_2^{p_2(\alpha_2 -1)+1}} 
\Big]^{1/p_2}
 = \Big[\frac{\Gamma(\frac{1}{4})}{2^{\frac{1}{4}}}  \Big]^{1/3}
 = \frac{\Gamma(\frac{1}{4})^{1/3}}{3^{1/12}}
\end{gather}
and since $m = 2$, $l = 2$, $d = m(l+1) + 2=8$, $\eta_{11} = \frac{3}{2}$, 
$\eta_{12} = \Gamma\big(1 + \frac{1}{q}\big) = \Gamma\big(\frac{7}{6}\big)$, 
$\alpha_1 = \frac{2}{3}$, $\alpha_2 = \frac{1}{2}$, we have
\begin{gather}\label{G1}
    G_1 = 8^5\cdot 8^6\Big (\frac{\Gamma(\frac{1}{3})^{1/2}}{2^{1/6}}
    \frac{1}{\Gamma(\frac{2}{3})}\Big)^{6}\Big(\frac{3}{2}\Big)^6
\Big(\frac{1}{3\cdot 2\cdot \frac{3}{2}}\Big)^\frac{3}{2}\frac{1}{6}, \\
\label{G2}
    G_2 =  8^5\cdot 8^6\Big( \frac{\Gamma(\frac{1}{4})^{1/3}}{3^{1/2}} 
\frac{1}{\Gamma(\frac{1}{2})}  \Big)^{12}\Gamma\Big(\frac{7}{6}\Big)^6
\Big(\frac{1}{3\cdot 2,2}   \Big )^{1/2}\frac{1}{2\cdot 6}, \\
D_2(z) = \frac{1}{1 - [C(z)D_1(z)]G_2} = \frac{1}{1 - d^qK^qz^qG_1G_2} 
= \frac{1}{1 - 8^6\cdot 8^6z^qG_1G_2},
\end{gather}
where $D_1(z) \equiv G_1$. Therefore
\[
\rho = \sup \big\{z \in \mathbb{R}: d^qK^q|z|^qG_1 < \frac{1}{G_2}\big\} 
= \sup \big\{z \in \mathbb{R}: 8^6\cdot 8^6|z|^6\cdot G_1 < \frac{1}{G_2}\big\},
\]
i.e.,
\[
\rho = \sup \big\{z \in \mathbb{R}: |z|
< \frac{1}{(8^6\cdot 8^6\cdot G_1\cdot G_2)^\frac{1}{6}}\big\},
\]
where $G_1, G_2$ are given by \eqref{G1} and \eqref{G2}, respectively, and
\[
 \Psi(w) = [C(w)D_1(w)D_2(w)]^{1/q} = \frac{dKw}{[1-d^qK^qw^qG_1G_2]^{1/q}}.
\]
By Theorem \ref{T},
\begin{equation}\label{SOE}
\|x(t)\| \leq \frac{dK\|x_0\|}{[1-d^qK^q\|x_0\|^qG_1G_2]^{1/q}}e^{-t} 
=  \frac{8\cdot 8\|x_0\|}{[1-8^6\cdot 8^6\|x_0\|^qG_1G_2]^{1/6}}e^{-t}
\end{equation}
for the solution $x(t)$ of the equation \eqref{ES} satisfying the
 initial condition $x(0) = x_0$ with $\|x_0\| < \rho$.

This means that for the solution $v(t)$ of the equation \eqref{EX}, which 
is the first coordinate of the solution $x(t) = (v(t), v'(t))$, 
the inequality \eqref{SOE} holds.


\subsection*{Acknowledgements}
The authors would like to express their gratitude to the anonymous referees 
for their valuable comments which significantly improved the original manuscript.

This work was supported by the Slovak Research and Development Agency 
under the contract APVV -14-0378 and by the Slovak Grant Agency VEGA-M\v S, 
project No. 1/0078/17.


\begin{thebibliography}{00}

\bibitem{BA} Bainov, D.; Simeonov, P.,
\emph{Integral Inequalities and Applications, Mathematics and Its Applications}
Vol. 57, Kluwer  Academic Press, Dortrecht, Boston, London 1992.

\bibitem{BT} Bagley, R. L.; Torvik, P.  J.;
\emph{On the appearance of the fractional derivative in the behavior of real materials},
J. Appl. Mech., \textbf{51} (1984), 294--298.

\bibitem{BRES} Brestovansk\'a, E.; Medve\v{d},  M.;
\emph{Exponential stability of solutions of a second order systems of 
integrodifferential equations with the Caputo-Fabrizio fractional derivatives},
Progr. Fract. Differ. Appl., \textbf{2(3)} (2016), 178--192.

\bibitem{DR} Delbosco, D.; Rodino, L.;
\emph{Existence and uniqueness for a nonlinear fractional differential equation},
J. Math. Anal. Appl., \textbf{204} (1996), 609--625.

\bibitem{FT} Furati, K. M.; Tatar, N.-E.;
\emph{Power type estimate for a nonlinear fractional differential equations},
Nonlinear Analysis, \textbf{62} (2015), 1025--1036.

\bibitem{FT1} Furati, K. M.; Tatar, N.-E.;
\emph{An existence result for a nonlocal fractional differential equations},
J. Fract. Calculus,  \textbf{26} (2004), 43--51.

\bibitem{KMT} Kirane, M., Medve\v{d}, M., Tatar, N.-E.;
\emph{Semilnear Volterra integrodifferential problems with fractional
 derivatives in the nonlinearities},
Abstr. Appl. Anal., Vol. \textbf{2011}, Art. ID 510314.

\bibitem{LAK} Lakshmikantham, V.; Leela, S.; Vasundhara, Devi J.;
\emph{Theory of Fractional Dynamical Systems},
Cambridge Scientific Publishers, Cambridge 2009.

\bibitem{LZ} Li, C. P. F.; Zhang, R.,
\emph{A survey on the stabiliy of fractional differential equations},
Eur. Phys. J., Special Topics \textbf{193} (2011), 27--47.

\bibitem{LD} Liu, L. L.; Duan ,J.~S.;
\emph{A detailed analysis for the fundamental solution of fractional 
vibration equation}, Open Math.,  \textbf{13} (2015), 826--838.

\bibitem{MAT1} Matigon, D.;
\emph{Stability result on fractional differential equations with applications 
to control processig}, Proceedings of IMACS-SMC, Vol. 2, Lille, France, 
July 1996, 963--968.

\bibitem{MAT2} Matigon, D.;
\emph{Stability properties for generalized fractional differential systems},
Proc. of Fractional Differential Systems: Models, Methods and Appl.,
 \textbf{5} (1998), 145--158.

\bibitem{LI} Li, C. P.; Zhang F.~ R.;
\emph{A survey on the stability of fractional differential equations},
The European Physical Journal, Special Topics, \textbf{193} (2011), 27--47.

\bibitem{MED1} Medve\v{d}, M.;
\emph{A new approach to an analysis of Henry type integral inequalities and their
Bihari type versions},
J. Math. Anal. Appl., \textbf{214} (1997), 349--366.

 \bibitem{MED2} Medve\v{d}, M.;
 \emph{Exponential stability of solutions of nonlinear differential equations 
with Riemann-Liouville fractional integrals in the nonlinearities},
 Proceedings of 4-th Scientific Colloquium, Prague June, 24--26 (2014), 10--20.

\bibitem{MED3} Medve\v{d}, M.;
\emph{Functional-differential equations with Riemann-Liouville integrals 
in the nonlinearities},
Mathematica Bohemica, \textbf{139(4)} (2014), 587--595.

\bibitem{MED4} Medve\v{d}, M.;
\emph{Singular integral inequalities with several nonlinearities and 
integral equations with singular kernels},
Nonlinear Oscillations, \textbf{11(1)} (2008), 71--80.

\bibitem{NA} Naber, M.,
\emph{Linear fractionally damped oscillator},
International Journal of Differential Equations,
\textbf{2010} (2010), Art. ID 197020, 12 pages.

\bibitem{PER} Perko, L.;
\emph{Differential Equations and Dynamical Systems},
Text in Applied Math. 7, Springer, New York 2001.

\bibitem{PET1} Petr\'a\v{s}, I.;
\emph{Stability of fractional-order systems with rational orders: Survay},
Fract. Calc. Appl. Anal., \textbf{12} (2009), 272--298.

\bibitem{PET2} Petr\'a\v{s}, I.;
\emph{Fractional-order Nonlinear Systems. Modeling, Analysis and Simulations},
Higher Education Press, Springer 2011.

\bibitem{PIN} Pinto, M.;
\emph{Integral inequality of the Bihari type and applications},
Funkcialaj Ekvacioj, \textbf{33} (1990), 387--401.

\bibitem{POD} Podlubn\'y, I.;
\emph{Fractional Differential Equations},
Academic Press, San Diego 1999.

\bibitem{PUT} Putzer, E. J.;
\emph{Avoiding the Jordan canonical form in the discussion of linear 
systems with constant coefficients}, Am. Math. Monthly, \textbf{73} (1966), 2--7.

\bibitem{SH} Seredynska, M.; Hanyga, A.;
\emph{Nonlinear differential equations with fractional damping with applications 
to the 1dof and 2dof pendulum}, Acta Mech., \textbf{176} (2005), 169-183.

\bibitem{TA} Tatar, N.-E.;
 \emph{Exponential decay for a system of equations with distributed delays},
 J.  Applied Math., \textbf{2015}, Art. ID 981383.

\bibitem{TA1} Tatar, N.-E.;
\emph{The existence of mild and classical solutions for a second-order 
abstract fractional problem}, Nonlinear Analysis, \textbf{73} (2010), 3130--3139.

\bibitem{TA2} Tatar, N.-E.;
\emph{Existence results for an evolution equation problem with fractional nonlocal conditions},
Computer and Mathemaics with Applications, \textbf{60} (2010), 2971--2982.

\bibitem{ZHO} Zhou, Y.;
\emph{Basic Theory of Fractional Differential Equations, World Scientific},
New Jersey, London Singapur, 2014.

\bibitem{ZHO1} Zhou, Y.;
\emph{Existence and uniqueness of solutions for a system of fractional
 differential equations}, Fractional Calculus \& Applied Analysis,
 \textbf{12(2)} (2009), 195--204.

\end{thebibliography}

\end{document}