\documentclass[12pt,twoside,a4paper]{article}
\usepackage[english]{babel}
\usepackage{pdfsync}
\newcommand{\thepages}{18--27}
\newcommand{\doinumber}{002}

\usepackage{umj_10_1}            % Style package of Ural Mathematical Journal

\usepackage{amsmath}            % If the package features are used
\usepackage{amsfonts,amssymb}   % If there is need to use math alphabets
                                % \mathbb{} - amsfonts.sty
                                % \mathfrak{} - eufrak.sty, amsfonts.sty
                                % symbols of font Euler - amssymb.sty
\newcommand{\Leftclt}{Abdulla Kh.~Abdullayev, Abdulla A.~Azamov,  Marks B.~Ruziboev}   % Left page header - name & family name of article author
\newcommand{\Rightclt}{An Explicit Estimate for Approximate Solutions} %% Right page header - Short title of article
\makeatletter\def\Hy@xspace@end{}\makeatother


\usepackage{enumitem}

%\usepackage{refcheck}
\setcounter{page}{18} %íóìåðàöèÿ ñ ...
\sloppy
\usepackage{srcltx}             % For easy move between WiNEdT and YAP

\begin{document}

\date{\today}

\Mainclt            % Page header print on the first page
% Enter article title and authors
% Reference to grants is a footnote to article title. You may use also Acknowledgments section for this purpose if you wish.

\begin{Titul}
{\large \bf \uppercase{AN EXPLICIT ESTIMATE FOR APPROXIMATE SOLUTIONS OF ODES BASED ON THE TAYLOR FORMULA}}%\footnote{This work was supported by ... (project No. NN-N). You may use also Acknowledgments section for this purpose if you wish.}
\\[3ex]
{{\bf Abdulla  Kh. Abdullayev} \\[2ex] {\small  V.I.~Romanovskiy Institute of Mathematics of Uzbekistan Academy of Sciences,\\ 9 University Str., 100174 Tashkent, Uzbekistan\\[5pt] \hrefemail{abduganiax@mail.ru} }\\[3ex]
{\bf Abdulla A. Azamov} \\[2ex] {\small V.I.~Romanovskiy Institute of Mathematics of Uzbekistan Academy of Sciences,\\ 9 University Str., 100174 Tashkent,   Uzbekistan \\[5pt] \hrefemail{abdulla.azamov@gmail.com} }\\[3ex]
{\bf Marks B. Ruziboev} \\[2ex] {\small Kimyo International University,\\ 156 Shota Rustaveli Str., 100121 Tashkent, Uzbekistan\\[5pt] \hrefemail{marx.ruziboev@gmail.com} }}\\[5ex]
\end{Titul}


\begin{Annot}
{\bf Abstract:}  In this paper, we consider a third-order explicit scheme based on Taylor's formula to obtain an approximate solution for the Cauchy problem of systems of ODEs. We prove an estimate for the accuracy of the approximate solution with an explicit constant that depends only on the right-hand side of the equation and the domain of the solution.

{\bf Keywords:} Dynamical systems, Cauchy problem, Approximate solution, Taylor formula, Accuracy of approximate solution, Level of accuracy, Error term.
\end{Annot}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %Please replace what follows by the body %%
 %of your article                         %%
                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction} %This section must have a number.
\setcounter{equation}{0}
It is needless to note the importance of estimating accuracy for approximate solutions of ODEs. Here we consider the problem
\begin{equation}\label{eq:ivp}
 \dot{x}(t)=f(x(t)),\quad x(0) = x_{0},
\end{equation}
where $x\in D \subset\mathbb{R}^{d}$ and $D$ is a convex domain. In what follows, we assume that the function $f: D\to  \mathbb R^{d}$  is three times differentiable with continuous derivatives in $D$. In practice, establishing the highest possible accuracy of an approximate solution is one of the key problems. Thus, the efficiency of an approximate solution is determined by its accuracy. Let $x_{*}(t)$ be a solution of \eqref{eq:ivp} in the interval $0\leq {t}\leq{T}$ for some $T>0$, and let $\widehat{x}(t)$  be its approximate solution (obtained by some scheme) on the same interval.  The accuracy of the scheme is expressed by an inequality of the form
\begin{equation}\label{eq:accdef}
\sup\limits_{0\leq t\leq T}\vert{x_{*}(t)-\widehat{x} (t)}\vert\leq C e^{LT}h^{s},
\end{equation}
where $L$ is the  Lipschitz constant of the function $f$, $h$ is the mesh size, and $s$ is the order of accuracy of the method. Approximate solution schemes can be implicit and explicit. In this paper, we are concerned with explicit methods. Frequently used one-step approximate solution schemes can be divided into two groups:
\vspace{-5pt}\begin{itemize}[label=\textup{(\alph*)}]\itemsep=-5pt
\item[(1)]	Schemes based on Taylor's expansion of the solution;
\item[(2)]	Runge--Kutta-type methods.
\end{itemize}
Schemes based on Taylor's expansion are easy to implement, but experts prefer Runge--Kutta-type methods. This preference is caused by the fact that the error estimates for schemes based on Taylor's formula contain derivatives of the function $f$, which can be challenging to estimate. However, with the development of computer algebra, computations of derivatives of a rather wide range of functions can be automated \cite{BuCoLo, GaGe, GeCzLa}. Therefore, schemes based on Taylor's formula can be implemented without any extra hurdles, due to the simplicity of implementation.

On the other hand, since the 1960s, computer simulations have been used extensively to study dynamical systems described by nonlinear systems of ordinary differential equations (ODEs). Recently, with the rise of computational power and computers being widely available, computer-assisted proofs have come into play. These proofs, however, require verification of the accuracy of the scheme, i.e., proof of inequality \eqref{eq:accdef} with explicit constants.

Estimates for the accuracy of approximate solutions to ODEs are studied extensively in literature \cite{Bab,BeZh,BrCaPe}. For example, the monographs \cite{Lam, Atk, HaNoWa} contain estimates for Runge--Kutta methods. Part~III of the well-known monograph \cite{Bak} is devoted to approximate solutions of ODEs. Schemes based on Taylor's formula are briefly discussed in the first section from a methodological viewpoint. But no estimates are provided. It is surprising that, among the vast amount of literature devoted to the approximate solutions of ODEs, we did not find estimates with explicit constants. Most estimates give the order of approximation, which is insufficient if we want to use approximate solutions in the proofs. Except \cite{HaNoWa}, which, citing \cite{Bieb}, gives an inequality for the Runge--Kutta method. In~\cite[Sect.~2, Part~II]{KrBoMo}, the authors mention a scheme based on Taylor's formula, but do not consider the problem of accuracy estimates. Some authors claim that if one takes the first $n$ terms of Taylor's expansion, then the error term will be of the form
$$\max\limits_{0\leq n\leq N}|x_{n}-y(x_{n})|=O(h^{n}).$$
But they neither provide proof nor speak about constants involved in $O(h^n)$.

In \cite{BuCoLo}, which is one of the most comprehensive monographs on approximate solutions of ODEs, the authors claim that the difference between exact and approximate solutions is estimated by the remainder term of Taylor's expansion and in just one step it will be $O(h^{p+1})$ \cite[Sect.~318, p.~180]{But}; no further details are given. In~\cite{Milbook, Milpap}, the problem of estimating the error is investigated for methods of approximation of the integral
$$
x(t+h)=x(t)+\int\limits_{t}^{t+h}f[x(t)]dt.
$$
Paper  \cite{But} provides explicit estimates for the approximate computation of this integral. This is equivalent to considering the first term of Taylor's expansion, which provides the first-order Taylor approximation scheme. The author considers multi-step approximate solution schemes but does not give inequalities of the form \eqref{eq:accdef}.

It turns out that schemes based on Taylor's formula are more convenient than schemes based on Runge--Kutta methods for obtaining explicit constants in \eqref{eq:accdef}. For example, in  \cite{AzAbTi}, an estimate of type \eqref{eq:accdef} is obtained for a second-order scheme based on Taylor's formula. Also, they are easier to implement in numerical approximations than Runge--Kutta-type methods. Keeping in mind an application of the approximate solution schemes in computer-aided proofs, in the present paper, we consider a third-order scheme based on Taylor's formula for Cauchy problem \eqref{eq:ivp} and give an explicit constant $C$, for which inequality \eqref{eq:accdef} holds with  $s=3$.

\section{Approximate solution schemes based on Taylor's formula}
\setcounter{equation}{0}
For a continuously differentiable function $f(x)$, the initial value problem \eqref{eq:ivp} has a unique solution; however, it is challenging to estimate the interval where the solution exists  \cite{Hart}. On the other hand, for approximation schemes, we need to know the existence of the solution. Therefore, our first standing assumption is the following.

\medskip
\textbf{Assumption A.} \emph{Fix $T>0$. The solution $x_{*}(t)$ to the Cauchy problem \eqref{eq:ivp} is defined on the interval $[0,T]$.
}

\medskip
Usually, when considering approximate solutions, one tries to find values of the approximate solution $x_{*} (t)$ on a mesh ${0 = t_0<t_1<t_2< \cdot\cdot\cdot <t_n=T}$. Here, to simplify the exposition, without loss of generality, we consider the uniform mesh $t_n=nh$, where $h={T}/{N}$, $N\in\mathbb N^{+}$.

We start with the Taylor expansion of the exact solution with accuracy $O(h^3)$:
\begin{equation*}
x_{*}(t+nh)=x_{*}(t)+h\dot{x}_{*}(t)+\frac{h^2}{2}\ddot{x}_{*}(t)+\frac{h^3}{6}\dddot{x}_{*}(t)+R_{4}(t,h).
\end{equation*}

By definition, we have
\begin{gather*}
\dot{x}_{*}(t)=f[{x}_{*}(t)],\\
\ddot{x}_{*}(t)=f'[{x}_{*}(t)]f[{x}_{*}(t)],\\
\dddot{x}_{*}(t)=f''[{x}_{*}(t)]f[{x}_{*}(t)]f[{x}_{*}(t)]+f'[{x}_{*}(t)]f'[{x}_{*}(t)]f[{x}_{*}(t)].
\end{gather*}
We also need the fourth derivative  $x^{IV} (t)$. To shorten the notation, we interpret derivatives of the function  $f$  as operators acting on $x_*(t)$ and write
\begin{gather*}
\dot{x}_{*}(t)=f[{x}_{*}(t)],\\
\ddot{x}_{*}(t)=(f'f)[{x}_{*}(t)],\\
\dddot{x}_{*}(t)=(f''ff+f'f'f)[{x}_{*}(t)].
\end{gather*}
In particular, on the uniform mesh,
\begin{gather*}
      x_{*}[(n+1)h]=x_{*}(nh)+hf[x_{*}(nh)]+\frac{h^2}{2}(f'f)[{x}_{*}(nh)] \\
      {}+\frac{h^3}{6}(f''ff+f'f'f)[{x}_{*}(nh)]+R_{4}(nh,h),\quad n=0,1,2,\ldots,N-1.
\end{gather*}
By neglecting the remainder term, we obtain the following recurrent formula for the approximate solution:
\begin{equation}\label{eq:asli}
 x_{n+1}=hf(x_{n})+\frac{h^2}{2}(f'f)(x_n)+\frac{h^3}{6}(f''ff+f'f'f)(x_n).
 \end{equation}

It is expected that elements of the sequence $x_n$ $(n=1,2,\ldots,N)$  defined by \eqref{eq:asli} are close to the values of the exact solution at the points  $x_* (h), x_* (2h),\ldots, x_* (Nh)$. An intuitive way to measure this closeness is to compute the value
 \begin{equation}\label{eq:maxdiff}
 \rho = \max_{1\leq n\leq N}|x_*(nh)-x_n|.
 \end{equation}
However, the quantity \eqref{eq:maxdiff} does not provide any information about the behavior of the solution on the interval $((n-1)h,nh).$  The main aim of this paper is to derive such estimates.  If we want to apply numerical solutions in computer-aided proofs, then we cannot ignore the behavior of the system on the interval $(nh,(n+1)h).$ Certainly, for  $t \in (nh,(n+1)h)$, we can estimate $x_* (t)$ as follows. For $t \in [0,T]$, let $n(t)=[{t}/{h}]$. Then $|t-n(t)h|<h$  and
 $$
 |x_*(t)-x_*(n(t)h)|\leq \int\limits_{n(t)h}^{h}|f[x(s)]|ds\leq M_0h,
 $$
 where $M_0$ is the maximum value of $|f(x)|$ on some compact subset of~$D$. The above inequality together with \eqref{eq:maxdiff} imply that
 $$
 |x_* (t)-x_n |\leq M_0 h+|x_* [n(t)h]-x_n |\le  M_0 h +\rho,
 $$
which  shows that the difference between the exact and approximate solutions is, at best, of order~$h$. However, this estimate is very rough and insufficient for our purposes. To get a better estimate, we use generalized Euler polygons.

 From now on, we fix $N,h,$ and $T>0$ such that $Nh = T$ and define
 $$
 \sigma_t = nh\quad \text{if}\quad t\in[nh, (n+1)h)\quad \text{for}\quad n = 0 , \ldots, N-1.
 $$
 We start with the definition of an approximate solution.
\begin{dfnen}A continuous function $\widehat{x} (t): \mathbb R \to [0:T]$ satisfying the equation
 \begin{equation}\label{eq:appsol}
 \widehat{x}(t)=x_0+\int\limits_{0}^{t}\Big(f[\widehat{x} (\sigma_s)]+(s-\sigma_s)(f'f)[\widehat{x} (\sigma_s)]+\frac{(s-\sigma_s)^2}{2}(f''ff+f'f'f)[\widehat{x} (\sigma_s)]\Big)ds
 \end{equation}
 is called an approximate solution of \eqref{eq:ivp}.
 \end{dfnen}

Although \eqref{eq:appsol} looks like an integral equation, it is a recurrent formula, and we can construct $\widehat{x}(t)$ explicitly step by step.  Therefore, \eqref{eq:appsol} defines a function  $\widehat{x}(t)$ on $[0, T]$ as an approximate solution; i.e., our notion of approximate solution is well defined.

\begin{lemen}
The equality $\widehat{x}(nh)=x_n$ holds for all $n > 0$.
\end{lemen}

The lemma is proved easily by induction on the intervals  $[0,nh]$.
We use this  lemma to reduce the problem to deriving estimates of the form \eqref{eq:accdef} for the difference $|x_* (t)-\widehat{x}(t)|.$

The remainder term of the Taylor expansion is given by the formula (see for example,  \cite{Car})
$$R_4(t,h)=\int\limits_{0}^{1}\frac{(1-s)^3}{3!}x_*^{(IV)}(t+sh)h^4ds,$$
where
$$x_*^{(IV)}(s)=\big(f'''fff+3f''f'ff+f'f''ff+f'f'f'f\big)[x_*(t)].$$
By the formula for $R_4$, the error term is estimated by the maximum values of the derivatives of $f$. To estimate the derivatives effectively along the solution and approximate solution, we need them to stay in some compact set. Therefore, we assume the following.

\medskip
\textbf{Assumption B.} \emph{Let $K$ be a convex and compact domain in  $\mathbb R^d$.  We assume that values of the exact solution $x_* (t)$ and the approximate solution $\widehat{x} (t)$ remain in $K$  for all $t\in [0, T].$}

\medskip
Since, $K$ is compact by definition,  Assumptions~$A$ and~$B$  allow us to define
\begin{gather}\label{eq:bounds}
\begin{gathered}
M_0=\max_{x\in{K}}|f(x)|,\quad M_1=\max_{x\in{K}}{\|f'(x)\|},\\
M_2=\max_{x\in{K}}{\|f''(x)\|},\quad  M_3 =\max_{x\in{K}}{\|f'''(x)\|}.
\end{gathered}
\end{gather}
Using these quantities, we obtain the following estimate for the remainder term of the Taylor formula:
$$
|R_4(nh,h)|\le\frac{h^4}{24}\,\max_{0\leq t\leq T}\,|x_*^{IV}(s)|ds\leq \frac{h^4}{24}\left(M_3M_0^3+4M_2M_1M_0^2+M_1^3M_0\right).
$$
In the literature, when estimating the accuracy of approximate solutions, many authors claim that the error is bounded by $|R_4(nh,h)|$. However,  this is not true,  since $R_4 (nh,h)$ is the difference between  $x_* (t)$  and its Taylor expansion, which does not directly imply any conclusions for the difference  $x_* (t)-\widehat x(t)$.
In the present paper, we show that it is possible to obtain an explicit estimate for the latter.
The main result of the this paper is the following theorem.

\begin{theoen}\label{thm}
Under Assumptions~$A$ and~$B,$ the following inequality holds\textup:
$$
|x_* (t)-\widehat x(t)|\leq \frac{e^{M_1T}-1}{6M_1}(L_0+L_1h+ L_2h^2)h^3,
$$
where
\begin{gather*}L_0=5M_0^2M_1M_2+M_0M_1^3+M_0^3M_3,\\
L_1=\frac{1}{4}(M_0^3M_2^2+4M_0^3M_1M_3+9M_0^2M_1^2M_2),\\
L_2=\frac{1}{2}(M_0^4M_2M_3+M_0^3M_1^2M_3+2M_0^3M_1M_2^2+2M_0^2M_1^3M_2).
\end{gather*}
\end{theoen}
Note that our estimate for the accuracy of the method is explicit and can be computed effectively in terms of the right-hand side of the initial value problem~\eqref{eq:ivp}.

\section{Proof of Theorem~\ref{thm}}
\setcounter{equation}{0}
In this section, we prove the main theorem. The key ingredient of the proof is a discretization of the time $t$ using a piecewise constant function $\sigma_s$.

In what follows, we repeatedly use the following formula for the derivative of the approximate solution. By \eqref{eq:appsol}, for $t\neq nh$ $(n=1,2,\ldots,N-1),$
\begin{equation}\label{eq:xshap}
\dot{\widehat x}(t)=f[\widehat x(\sigma_t)]+(t-\sigma_t)(f'f)[\widehat x(\sigma_t)]+\frac{(t-\sigma_t)^2}{2}(f''f^2+(f')^2f)[\widehat x(\sigma_t)].
\end{equation}
Further, taking into account that  $0\leq t-\sigma_t \leq h$, we have
\begin{equation}\label{eq:|xshap|}
|\dot{\widehat x}(\sigma_t)|\leq M_0 \Big(1+hM_1+\frac{h^2}{2}(M_2M_0+M_1^2)\Big).
\end{equation}


\proofen For the exact solution, we have the equality
$$x_*(t)=x_0+\int\limits_{0}^{t}f[x_*(s)]ds.$$

Using this equality and equation \eqref{eq:appsol} and adding and subtracting the term $$\int_0^tf[\widehat x(s)]ds,$$ we obtain
$$
|x_* (t)-\widehat x(t)|\leq \int\limits_{0}^{t}I(s)ds+ \int\limits_{0}^{t}\big|f[x_*(s)]-f[\widehat x(s)]\big|ds,
$$
where
$$I(s)=f[\widehat x(s)]-f[\widehat x(\sigma_s)]-(s-\sigma_s)(f'f)[\widehat x(\sigma_s)]-\frac{(s-\sigma_s)^2}{2}(f''f^2+f'^{2}f)[\widehat x(\sigma_s)].$$

We are going to derive an upper estimate of the form $C h^3$ for $I(s)$. By the fundamental rule of the calculus, we have
$$f[\widehat x(s)]-f[\widehat x(\sigma_s)]=\int\limits_{\sigma_s}^{s}\frac{df[\widehat x(r)]}{dr}dr=\int\limits_{\sigma_s}^{s}f'[\widehat x(r)]\dot{\widehat x}(r)dr.$$
Substituting the derivative of the approximate solution $\dot{\widehat x}(r)$ given in  \eqref{eq:|xshap|} into the right-hand side, we obtain
 \begin{gather*}
 f[\widehat x(s)]{-}f[\widehat x(\sigma_s)]{=}\int\limits_{\sigma_s}^{s}f'[\widehat x(r)]\Big\{f[\widehat x(\sigma_r)]{+}(r{-}\sigma_r)(f'f)[\widehat x(\sigma_r)]
 +\frac{(r-\sigma_r)^2}{2}(f''f^2{+}f'^{2}f)[\widehat x(\sigma_r)] \Big\}dr.
 \end{gather*}

 Note that $\sigma_r =\sigma_s$ in this equation since $r\in[\sigma_s, s]$ and $\sigma$ is piecewise constant. Denote the latter term on the right-hand side by $C_1$:
 $$
 C_1=\int\limits_{\sigma_s}^{s}\frac{(r-\sigma_s)^2}{2}(f''f^2+f'^{2}f)[\widehat x(\sigma_r)]f'[\widehat x(r)]dr.
 $$
Using \eqref{eq:bounds}, we obtain the estimate
\begin{equation}\label{eq:estc1}
|C_1|\leq (M_2M_1M_0^2+M_1^3M_0)\int\limits_{\sigma_s}^{s}\frac{(r
-\sigma_s)^2}{2}dr\leq D_1h^3,
\end{equation}
where $D_1=(M_2M_1M_0^2+M_1^3M_0)/6.$

Therefore,  we obtain
\begin{gather}\label{eq:defi}
\begin{gathered}
I(s)=\int\limits_{\sigma_s}^{s}\big\{f'[\widehat x(r)]-f'[\widehat x(\sigma_s)]\big\}f[\widehat x(\sigma_s))]dr\\
+\int\limits_{\sigma_s}^{s}f'[\widehat x(r)](f'f)[\widehat x(\sigma_s)]dr -\frac{(s-\sigma_s)^2}{2}(f''f^2+f'^{2}f)[\widehat x(\sigma_s)]+C_1.
\end{gathered}
\end{gather}
Denote the first term of this expression by $J(s)$:
$$ J(s)=\int\limits_{\sigma_s}^{s}\big\{f'[\widehat x(r)]-f'[\widehat x(\sigma_s)]\big\}f[\widehat x(\sigma_s)]dr.$$
We estimate $J(s)$, using the fundamental rule of calculus and the derivative of the approximate solution given by \eqref{eq:|xshap|}:
\begin{gather}\label{eq:defj}
\begin{gathered}
J(s)=\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}\frac{df'[\widehat x(u)]}{du}duf[\widehat x(\sigma_s)]dr=
\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}f''[\widehat x(u)]\Big\{ f[\widehat x(\sigma_u)]
\\
{}+(u-\sigma_u)(f'f)[\widehat x(\sigma_u)]+\frac{(u-\sigma_u)^2}{2}(f''f^2+(f')^2f)[\widehat x(\sigma_u)]\Big\}f[\widehat x(\sigma_s)]dudr.
\end{gathered}
\end{gather}
Define
\begin{equation}\label{eq:defc2}
\begin{split}
C_2=\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}\Big\{ (u-\sigma_u)f''[\widehat x(u)](f'f)[\widehat x(\sigma_u)]f[\widehat x(\sigma_s)]\\
{}+\frac{(u-\sigma_u)^2}{2}f''[\widehat x(u)](f''f^2+(f')^2f)[\widehat x(\sigma_u)]f[\widehat x(\sigma_s)]\Big\}dudr.
\end{split}
\end{equation}
Thus, taking into account that $\sigma_u=\sigma_s$ and using \eqref{eq:bounds}, we obtain the inequality
\begin{equation}\label{eq:estc2}
|C_2|\leq\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}\left|M_2M_1M_0^2(u-\sigma_s)+(M_2^2M_0^3+M_2M_1^2M_0^2)\frac{(u-\sigma_s)^2}{2}\right|dudr
\leq D_2h^3+D_3h^4,
\end{equation}
where
$$D_2=\frac{1}{6}M_1M_2M_0^2,\quad  D_3=\frac{1}{24}(M_2^2M_0^3+M_2M_1^2M_0^2).$$
Consequently, substituting \eqref{eq:defc2} into \eqref{eq:defj} and then  \eqref{eq:defj} into \eqref{eq:defi} and denoting $C_1+C_2$ by $C_3$, we obtain
\begin{gather}\label{eq:finali}
\begin{gathered}
I(s)=\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}f''[\widehat x(u)]f[\widehat x(u)]f[\widehat x(\sigma_s)]dudr+\int\limits_{\sigma_s}^{s}f'[\widehat x(r)](r-\sigma_r)(f'f)[\widehat x(\sigma_r)]dr\\
{}-\frac{(s-\sigma_s)^2}{2}(f''f^2+f'^{2}f)[\widehat x(\sigma_s)]+C_3=A(s)+B(s)+C_3,
\end{gathered}
\end{gather}
where we used the following notation:
\begin{gather*}
A(s)= \int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}\big\{ (f''f)[\widehat x(u)]-(f'f)[\widehat x(\sigma_s)]\big\}f[\widehat x(\sigma_s)]dudr,\\
B(s)=\int\limits_{\sigma_s}^{s}(r-\sigma_s)\big\{ f'[\widehat x(r)]-f'[\widehat x(\sigma_s)]\big\}(f'f)[\widehat x(\sigma_s)]dr.
\end{gather*}
Combining \eqref{eq:estc1} and \eqref{eq:estc2},  we obtain
\begin{equation}\label{eq:estc3}
| C_3 |\leq D_1 h^3+D_2 h^3+D_3 h^4.
\end{equation}
It remains to estimate $A(s)$ and $B(s)$.
We have
$$ A(s)=\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}\bigg(\int\limits_{\sigma_s}^{u} \frac{d}{dv}\left(f''[\widehat x(v)]f[\widehat x(v)]\right)dv\bigg)f[\widehat x(\sigma_s)]dudr.
$$
which implies
\begin{equation}\label{eq:estas}
A(s)=\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}\bigg(\int\limits_{\sigma_s}^{u}\big\{ f'''[\widehat x(v)]\dot{\widehat x}(v)f[\widehat x(v)]+f''[\widehat x(v)]f'[\widehat x(v)]\dot{\widehat x}(v)\big\}dv\bigg)f[\widehat x(\sigma_s)]dudr.
\end{equation}

Therefore, taking the absolute value of the expression under the outer integral (integration with respect to $r$), using estimates \eqref{eq:bounds} and \eqref{eq:|xshap|}, and taking into account that $t-\sigma_t\le h$ and $\sigma_u=\sigma_s = \sigma_r$, we obtain
%$$ |A(s)|\leq \int\limits_{\sigma_s}^{s}\left| \int\limits_{\sigma_s}^{r}\left( \int\limits_{\sigma_s}^{u}\{f^{'''}[\widehat x(v)]\widehat x^{'}(v)f[\widehat x(v)]+f^{''}[\widehat x(v)]f^{'}[\widehat x(v)]\widehat x^{'}(v)\}dv\right)f[\widehat x(\sigma_s)]du \right|dr \leq$$
\begin{gather*}\label{eq:as}
\begin{gathered}
|A(s)|
\le\int\limits_{\sigma_s}^{s}\int\limits_{\sigma_s}^{r}\int\limits_{\sigma_s}^{u}\Big(M_0+hM_1M_0+\frac{h^2}{2}(M_2M_0^2+M_1^2M_0)\Big)
(M_3M_0^2+M_0M_1M_2)dvdudr \\
 \leq(M_3M_0^2+M_0M_1M_2)\Big(M_0+hM_1M_0+\frac{h^2}{2}(M_2M_0^2+M_1^2M_0)\Big)\frac {h^3}{6}.
\end{gathered}
\end{gather*}
Similarly,
\begin{gather*}
B(s)=\int\limits_{\sigma_s}^{s}(r-\sigma_s)\bigg\{\int\limits_{\sigma_s}^{r}\frac{d}{du}f'[\widehat x(u)]du \bigg\}(f'f)[\widehat x(\sigma_s)]dr\\
=\int\limits_{\sigma_s}^{s}(r-\sigma_s)\int\limits_{\sigma_s}^{r}f''[\widehat x(u)]\dot{\widehat x}(u)du \}(f'f)[\widehat x(\sigma_s)]dr.
\end{gather*}
Again, using \eqref{eq:bounds} and \eqref{eq:|xshap|}, and taking into account that $t-\sigma_t\le h$ and $\sigma_s = \sigma_r$ in the above equation, we obtain
% $$ |B(s)|\leq \int\limits_{\sigma_s}^{s}\left|(r-\sigma_s)\int\limits_{\sigma_s}^{r}f^{''}[\widehat x(u)]\widehat x^{'}(u)du(f^{'}f)[\widehat x(\sigma_s)]\right|dr\leq$$
\begin{gather}\label{eq:estbs}
\begin{gathered}
|B(s)|\leq \frac{(r-\sigma_s)^3}{3}M_0M_1M_2\Big(M_0+hM_1M_0+\frac{h^2}{2}(M_2M_0^2+M_1^2M_0)\Big)\\
\leq M_0M_1M_2\Big(M_0+hM_1M_0+\frac{h^2}{2}(M_2M_0^2+M_1^2M_0)\Big)\frac{h^3}{3}.
\end{gathered}
\end{gather}
Finally, substituting   \eqref{eq:estc3},  \eqref{eq:estas}, and \eqref{eq:estbs} into \eqref{eq:finali}, we obtain the inequality
\begin{gather*}
 |I(s)|\leq\frac{1}{6}(M_2M_1M_0^2+M_1^3M_0))h^3+\frac{1}{6}M_1M_2M_0^2h^3+\frac{1}{24}(M_2^2M_0^3+M_2M_1^2M_0^2)h^4 \\
 +(M_3M_0^2+M_0M_1M_2)\left(M_0+hM_1M_0+\frac{h^2}{2}(M_2M_0^2+M_1^2M_0)\right)\frac {h^3}{6}\\
 +M_0M_1M_2\left(M_0+hM_1M_0+\frac{h^2}{2}(M_2M_0^2+M_1^2M_0)\right)\frac{h^3}{3}\\
=\frac{1}{6}(L_0+L_1h+L_2h^2)h^3,
\end{gather*}
with the required constants $L_0$, $L_1$, and $L_2$.

Now, we %apply Hadamard lemma to  the difference $f[x_* (s)]-f[\widehat x(\sigma_s)]$.
use the compactness of the domain $K$ and smoothness of $f$ to obtain
\begin{equation}\label{eq:lipf}
\begin{split}
      |f[x_* (s)]-f[\widehat x(\sigma_s)]|
     %=\left|\int\limits_{0}^{1}\frac{\delta[rx_*(s)+(1-r)\widehat x(s)]}{\delta x}dr\{ x_*(s)-\widehat x(s)\}\right|\leq \\
       \le M_1|x_*(s)-\widehat x(s)|.
\end{split}
\end{equation}

Using inequalities \eqref{eq:xshap}, \eqref{eq:|xshap|}, and \eqref{eq:lipf}, we get the inequality
\begin{gather}\label{eq:diff}
\begin{gathered}
      |x_* (t) -\widehat x(t)|\leq \int\limits_{0}^{t}|I(s)|ds+\int\limits_{0}^{t}|f[x_* (s)]-f[\widehat x(\sigma_s)]|ds \\
      \leq\frac{1}{6}(L_0+L_1h+L_2h^2)h^3t+\int\limits_{0}^{t}M_1|x_*(s)-\widehat x(s)|.
\end{gathered}
\end{gather}
In  \eqref{eq:diff}, we apply Gr\"onwall's inequality. For our purposes, the following version is the most convenient. Let  $u:R \to R$ be a continuous function such that  $u(t)\ge 0$ for $t\ge 0$  and  $$u(t)\leq Ct+M \int\limits_{0}^{t}u(s)ds$$ for some $C, M>0$. Then, the following inequality holds:
$$u(t)\leq C\frac{e^{Mt}-1}{M}.$$
Applying Gr\"ownall's inequality to estimate  \eqref{eq:diff}, we obtain
$$ |x_* (t) -\widehat x(t)|\leq \frac{1}{6}(L_0+L_1h+L_2h^2)h^3\frac{e^{M_1t}-1}{M_1}.$$
This completes the proof. \hfill$\square$


\section{Conclusion}
\setcounter{equation}{0}
\begin{enumerate}
  \item  If a nonautonomous system is considered in a $d$-dimensional space, then we can interpret it as an autonomous system in the $(d+1)$-dimensional space.  In particular, we can consider the following Cauchy problem:
\begin{equation*}
\dot{x}=f(t,x),\quad  x(0)=x_0 \Leftrightarrow\left\{
\begin{array}{ll}
\dfrac{dx}{dt}=f(\xi,x), & x(0)=x_0,\\[1ex]
\dfrac{d\xi}{dt}=1,& \xi(0)=0.
\end{array}\right.
\end{equation*}
Therefore, the difference between the exact and approximate solutions can be estimated by the same expression with the constant $\sqrt{M_0^2+1}$ insted of~$M_0$.
\item
In the proof of the estimate for the difference between the exact and approximate solutions, we obtained
$$\max\limits_{0\leq t\leq T}|x_* (t) -\widehat x(t)|\leq C\frac{e^{M_1t}-1}{M_1}h^3 $$
with the coefficient that is considerably larger than expected, where the constant $C$ is a fifth-order polynomial of the constants $M_0$, $M_1$, $M_2$, and $M_3$.  If these constants are not very large, i.e., on the order of 1, then the coefficient $C$  does not affect the choice of the mesh size. In this case, the mesh size would mostly depend on  ${(e^{M_1t}-1)}/{M_1}$. On the other hand, if the constants $M_0$, $M_1$, $M_2$, and $M_3$ are on about $10$, then the coefficient at $h^3$ is on the order of $10^5$, which would affect the choice of $h$ essentially, in certain cases, it may even invalidate the approximate solution scheme.
\item Another interesting question is whether it is possible to simplify the proof of the main theorem. The authors think that the proofs cannot be simplified considerably.
\item It is possible to prove a similar theorem with the order of accuracy $h^4$; i.e.,  we can consider the fourth-order scheme
\begin{gather*}
x_{n+1}=x_n+hf(x_n)+\frac{h^2}{2}(f'f)(x_n)+\frac{h^3}{6}\big(f''ff+f'f'f\big)(x_n)+\\
+\frac{h^4}{24}\big(f'''fff+3f''f'ff+f'f''ff+f'f'f'f\big)(x_n)
\end{gather*}
and prove an analogous theorem.
\end{enumerate}

\begin{Biblioen}
\bibitem{AzAbTi} Azamov A.\,A., Abdullaev A.\,X., Tilavov A.\,M.  On the derivation of an inequality for estimating the accuracy of an approximate solution to the initial value problem. \emph{Bull. Inst. Math}., 2022. Vol.~5, No.~5. P.~105--111.
\bibitem{Bab} Babenko K.\,I. \emph{Foundations of Numerical Analysis} [Osnovy chislennogo analiza]. 2nd edition. Moscow, Izhevck: Regulyarnaya i haoticheskaya dinamika, 2002. 848~p. (in~Russian)
\bibitem{Bak} Bakhvalov N.\,S. \emph{Numerical Methods (Analysis, Algebra, Ordinary Differential Equations)} [CHislennye metody (analiz, algebra, obyknovennye differencial'nye uravneniya)]. Moscow: Nauka, 1975. 632~p.
\bibitem{BeZh} Berezin I.\,S., Zhidkov P.\,N. \emph{Computing Methods}.  Vol.~2. Pergamon Press, 1965. 464~p.
\bibitem{Bieb} Bieberbach L. On the remainder of the Runge--Kutta formula in the theory of ordinary differential equations. \emph{J. Appl. Math. Phys. (ZAMP)}, 1951. Vol.~2. P.~233--248.
   \hrefdoi{10.1007/BF02579687} (in German)
 \bibitem{BrCaPe} Brenan K.\,E., Campbell S.\,L., Petzold L.\,R. \emph{Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations.} Philadelphia: SIAM, 1995.  XII+251~p. \hrefdoi{10.1137/1.9781611971224}
\bibitem{BuCoLo} Buchberger B., Collins G.\,E., Loos R., Albrecht R. \emph{Computer Algebra: Symbolic and Algebraic  Computation}.  Vienna: Springer,  1983. VII+283~p. \hrefdoi{10.1007/978-3-7091-7551-4}
    %Cite this chapter
    %Loos, R. (1983). Computing in Algebraic Extensions. In: Buchberger, B., Collins, G.E., Loos, R., Albrecht, R. (eds) %Computer Algebra. Computing Supplementa, vol 4. Springer, Vienna. https://doi.org/10.1007/978-3-7091-7551-4_12

\bibitem{But} Butcher J.\,C. \emph{Numerical Methods for Ordinary Differential Equations}. 3rd edition. United Kingdom: John Wiley \& Sons Ltd., 2016. 544~p.
\bibitem{Car} Cartan H. \emph{Formes Diff\'{e}rentielles}. Hermann Paris, 1967. 190~p. (in~French)
\bibitem{HaNoWa} Hairer E., Norsett S., Wanner G. \emph{Solving Ordinary Differential Equations I. Nonstiff Problems.} 2nd edition.  Berlin, Heidelberg: Springer, 2008. XV+528~p. \hrefdoi{10.1007/978-3-540-78862-1}
\bibitem{Hart} Hartman P. \emph{Ordinary Differential Equation}. SIAM, 2002. XVIII+624~p. \hrefdoi{10.1137/1.9780898719222}
\bibitem{GaGe} Gathen J., Gerhard J. \emph{Modern Computer Algebra}. 3rd edition. Cambridge University Press, 2013. 796~p. \hrefdoi{10.1017/CBO9781139856065}
\bibitem{GeCzLa} Geddes K.\,O., Czapor S.\,R., Labahn G. \emph{Algorithms for Computer Algebra}. NY: Springer,   2007. XXII+586~p. \hrefdoi{10.1007/b102438}
\bibitem{Atk} Kendall A.\,E.  \emph{An Introduction to Numerical Analysis}. 2nd edition. John Wiley \& Sons, Inc., 1989. 693~p.
\bibitem{KrBoMo} Krilov V.\,I., Bobkov V.\,V., Monastirniy P.\,I. \emph{Computational Methods} [Vychislitel'nye metody]. Vol. II. Moscow: Nauka, 1977. 399~p. (in~Russian)
\bibitem{Lam} Lambert J.\,D. \emph{Numerical Methods for Ordinary Differential Systems: The Initial Value Problem}. Wiley, 1992. 304~p.
\bibitem{Milbook} Milne W.\,E. \emph{Numerical Calculus}. Princeton University Press,  1949. 404~p.
\bibitem{Milpap} Milne W.\,E.  The Remainder in linear methods of approximation. \emph{J. Res. Nat. Bur. Stand.}, 1949. Vol.~43. Art. no.~RP2401. P.~501--511.
\end{Biblioen}

\end{document}
