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

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

\begin{document}
\title[\hfilneg EJDE-2012/11\hfil Cauchy-Kowalevski and polynomial ODEs]
{Cauchy-Kowalevski and polynomial ordinary differential equations}

\author[R. J. Thelwell, P. G. Warne,  D. A. Warne\hfil EJDE-2012/11\hfilneg]
{Roger J. Thelwell, Paul G. Warne, Debra A. Warne}  % in alphabetical order

\address{Roger J. Thelwell \newline
Department of Mathematics and Statistics,
James Madison University, Harrisonburg, VA 22807, USA}
\email{thelwerj@jmu.edu}

\address{Paul G. Warne \newline
Department of Mathematics and Statistics,
James Madison University, Harrisonburg, VA 22807, USA}
\email{warnepg@jmu.edu}

\address{Debra A. Warne \newline
Department of Mathematics and Statistics,
James Madison University, Harrisonburg, VA 22807, USA}
\email{warneda@jmu.edu}

\thanks{Submitted October 11, 2010. Published January 17, 2012.}
\subjclass[2000]{34A12, 34A34, 35A10}
\keywords{Automatic differentiation; power series; Taylor series;
 \hfill\break\indent polynomial ODE; majorant; error bound}

\begin{abstract}
 The Cauchy-Kowalevski Theorem is the foremost result guaranteeing
 existence and uniqueness of local solutions for analytic
 quasilinear partial  differential equations with Cauchy initial data.
 The techniques of Cauchy-Kowalevski may also be applied to
 initial-value ordinary differential equations.
 These  techniques, when applied in the polynomial ordinary
 differential equation setting, lead one naturally to a method
 in which coefficients of the series solution are easily
 computed in a recursive manner, and an explicit majorization
 admits a clear a priori error bound.
 The error bound depends only on immediately
 observable quantities of the polynomial system; coefficients,
 initial conditions, and polynomial degree.
 The numerous benefits of the polynomial system are shown for
 a specific example.
\end{abstract}

\maketitle
\allowdisplaybreaks

\section{Introduction}

The Cauchy-Kowalevski Theorem is the main tool in showing the
existence and uniqueness of local solutions for analytic quasilinear
partial differential equations (PDE) with Cauchy initial data.
Cauchy developed a proof in a restricted setting by 1842 \cite{Cauchy},
and in 1875 Kowalevski presented the full result \cite{Kowalevski};
existence of a unique solution to the general quasilinear system of
partial differential equations given initial conditions prescribed on
some non-characteristic curve. In \cite{Folland}, a proof in the fully
nonlinear setting is presented.  The Cauchy-Kowalevski argument is based
on the construction of a power series solution, in which the coefficients
of the series expansion are reconstructed recursively, and the method
of majorants applied to verify that this solution converges locally.
Convergence is demonstrated by comparison with the analytic solution of
an associated PDE.

Although the Picard-Lindel\"of Theorem is the fundamental local
existence argument for a large class of initial value ordinary
differential equations (IVODE), in 1835 Cauchy demonstrated existence
and uniqueness in the ODE setting, applying a majorant based argument
similar to that both
he and Kowalevski would later use in the PDE setting.
That is, Cauchy methods can be used to show that $u$ satisfies
the real analytic ODE $d_t u(t) = f(u(t))$, where $u(0) = u_0$
using a constructive approach, provided $f(u)$ is analytic near $u_0$.
A nice treatment may be found in \cite{Driver}.

Given that the power series solution is directly accessible
via the Cauchy-Kowa\-levski construction but that the method is rarely
applied suggests practical difficulties.
In fact, the coefficients of the series solution can be tedious
to construct as typically posed,
as is a key constant in the comparison solution.
In this paper, we demonstrate that a subtle
recasting of the ODE system meliorates these difficulties:
the coefficients of the analytic solution become
remarkably easy to recover, and a computable choice of the
key constant in the majorization leads to an \textit{attractive a priori error bound.}
To make these ideas clear, we consider the simple  quasilinear problem
\begin{equation}
\label{E:IVODE}
d_t u(t) = f(u(t)) := \frac{\mathrm{e}^{2 u(t)^2}}{\sin u(t)} ,
\quad\text{with } u(0) = 1.
\end{equation}

We first consider \eqref{E:IVODE} using the methods of Cauchy,
and identify steps in which the construction of solution becomes tedious.
We then recast the problem as a polynomial system, as might be done
when using Taylor series based automatic differentiation, and apply the
same methods.  It will be clear the computations necessary to
generate the series solution are basic, and that a simple majorization
which depends only on the magnitude of the \textit{initial conditions}, the \textit{degree}
of the polynomial system and the magnitude of the \textit{constant coefficients}
of the system leads to an explicit bound of the remainder when approximating with the Taylor Polynomial.
Although not demonstrated here, the method applied is quite general.
The authors view this note as complimentary to \cite{Warne2006}.
Most importantly, we conjecture that it may be possible to extend the method to analytic IVPDE.

\section{Recasting (non)linear ODE(s) as polynomial systems: Why?}

Ordinary differential equations, particularly nonlinear and those with
singularities, play a fundamental role in understanding the principles
that govern the world around us.
Left in their original
(or \textit{classic}) form, various analytic and numeric methods exist which are
problem specific, yet there remains a need for a systematic method to calculate
solutions of general problems.
The approach presented here, perhaps first introduced by Cauchy and
subsequently rediscovered and coupled with power series
methods by Fehlberg in 1964 (\cite{Fehlberg1964}) and others since, is
simple and surprisingly general.  A recasting of the original
ODE as a system of constant coefficient polynomial ODEs
via an introduction of auxiliary variables leads to a straight
forward iterative calculation of power series coefficients.
This allows a clear and systematic construction of numeric solutions and
provides an immediate and explicit \textit{a priori} error bound.

A polynomial system is useful computationally, and methods are
available to recast an impressively wide variety of ODEs into an augmented
polynomial system.
 This includes ODEs with right-hand sides involving compositions
of exponential, logarithmic, and trigonometric functions,
as well as those involving the algebraic operations, including
exponentiation of complex power and encompassing general reciprocals and singularities.
For example, the second order ODE,
\begin{equation}
\label{deb}
y''\Big[ {1 + \frac{{\sqrt 2 }}
{{{{y'}^2}}}{{\big( {\frac{x}
{{yy'}}} \big)}^{\sqrt 2  - 1}}} \Big]
= \Big[ {\sqrt 2 {{\big( {\frac{x}
{{yy'}}} \big)}^{\sqrt 2 }} + 1} \Big]
\Big[ {\frac{y}{{{x^2}}} - \frac{{y'}}
{x}} \Big] + \frac{{{\pi ^2}}}{{16}}y ,
\end{equation}
with ${( {})' }: = \frac{d}{{dx}}$, used to
model the torsional deformation of a compressible elastic solid
cylinder composed of a generalized Blatz-Ko material, does not appear amiable to classic power
series methods, and yet its series solution to arbitrary order
is easily computable from its equivalent polynomial system \cite{Paullet}.
While such computational advantage is not the focus here,
see \cite{Barrio, Griewank,Jorba2005,Rall} for examples and discussion in
the automatic differentiation setting or
\cite{Parker1996,Pruett2003,Perlin1964,Stewart2009} in the ODE setting.
Recasting an ODE as a polynomial system also makes analysis tractable. We now explore
one such application: a proof of existence and uniqueness.

\section{Cauchy solution: the classic setting}

We begin with the precarious assumption that a locally analytic solution $u(t)$ to \eqref{E:IVODE} exists, and repeatedly differentiate the equation, using the fact that $f(u)$ is analytic in $u$ near the initial condition.
\begin{align*}
d_t^2 u(t) & = d_u f(u) d_t u    \\
   & = -{\frac {{{\rm e}^{4\,{u}^{2}}} \left( 4\,u\sin u -\cos
 u  \right) }{ \left( -1+ \cos^{2} u \right) \sin u }},
 \notag \\
d_t^3 u(t) & = d_u^2 f(u) [d_t u]^2 + d_u f(u) d_t^2 u \notag \\
   & = -{\frac {{{\rm e}^{6\,{u}
^{2}}} \left( 32\,{u}^{2}\cos^{2} u+2\, \cos^{2} u + 16\,u\sin u \cos u -5-32\,{u}^{2} \right) }{ \left( 1-2\, \cos^{2} u +
  \cos^{4} u  \right) \sin u }}
 \\
d_t^4 u(t) & = d_u^3 f(u) [d_t u ]^3 + 3 d_u^2 f(u) d_t^2 u d_t u + d_u f(u) d_t^3 u \notag \\
  & = {\frac {{{\rm e}^{8\,{u}^{2}}} \left( M_1 + M_2  \right)  }{ \left( -1+3\, \cos^{2} u-3\, \cos^{4} u+ \cos^{6} u \right) \sin u }}
 \end{align*}
where
\begin{gather*}
M_1  =  -288\,{u}^{2} \cos^{3} u-22\, \cos^{3} u+384\,{u}^{3}
\sin u  \cos^{2} u -140\,u\sin u  \\
M_2  =  40\,u\sin u \cos^{2} u+37\,\cos u +288\,\cos u {u}^{2}
+384\,{u}^{3}\sin u
\end{gather*}
and
\begin{equation} \label{E:polynomial}
 d_t^n u(t)  = p_n(f(u),d_u f(u),d_u^2 f(u), \ldots , d_u^{n-1} f(u)),
\end{equation}
where $p_n(\cdot)$ denotes a polynomial in $n$ variables
(here taken from the set of derivatives of
$f$ with respect to $u$ of order less than $n$; i.e.,
 $\{ d_u^{k-1} f \}, k=1,\ldots,n$, and having positive integer
coefficients).
By this process, all coefficients of the power series representation of
$u(t)$ may be built;
\begin{equation}
\label{E:PSS}
u(t) = \sum_{k=0}^\infty \frac{1}{k!} d_t^{k} u(0) \,  t^k.
\end{equation}
Note that the form of the polynomial $p_n$ in expression
\eqref{E:polynomial}
allows the coefficients of the power series to be recovered
recursively, although the complexity of calculation may
(and usually does) grow exponentially.

By its very construction, this power series \eqref{E:PSS}
yields a unique classical solution to the initial-value ODE if it can be shown
to converge. Cauchy demonstrated convergence by comparison with
a related analytic initial-value ODE,
whose individual coefficients majorize (absolutely bound) those of
\eqref{E:PSS}.  We briefly illustrate the argument.
We begin with the assumption of the theorem that $f(u)$ is analytic
in some interval of radius $R \in \mathbb{R}$ about $u=1$,
and remark that in practice $R$ might be quite difficult to determine.
Then for any positive $r<R$, there exists
\[
C_\infty := \max_k \{|C_k|\} < \infty, \quad\text{where }
C_n = \frac{1}{n!} d_u^n f(1) r^n,
\]
which provides the bound
\[
\max_k \big| \frac{1}{k!} d_u^k f(1) \big| \le C_\infty r^{-k}
\]
on the Taylor coefficients of $f(u)$ about $u(0)=1$.
Next we define $g$ via the geometric series
\[
g(v) := \sum_{k=0}^\infty C_\infty r^{-k} (v-1)^k
= C_\infty \frac{r}{r-(v-1)} \quad\text{when } |v-1|<r,
\]
and the comparison initial-value ODE
\begin{equation}
\label{E:MIVODE}
d_t v(t) = g(v(t)) \quad\text{with } v(0) = 1.
\end{equation}
The form of equation \eqref{E:MIVODE} is motivated by the observation
that the polynomial $p_n$ generated in this case is identical in
form to that of \eqref{E:polynomial}, allowing a direct comparison
of coefficients of $u(t)$ with those of $v(t)$.
Also, $g(v)$ majorizes $f(u)$ near 1 and allows \eqref{E:MIVODE}
an analytic solution $v(t)$ near 0.
When $|v-1|<r$,
\[
|d_u^n f(1)| = n! \big| \frac{1}{n!} d_u^n f(1) \big|
\le n!\, C_\infty r^{-n} = d_v^n g(1)
\]
for all $n$.  Noting that the structure of the polynomial
in \eqref{E:polynomial}
is identical in \eqref{E:IVODE} and \eqref{E:MIVODE},
it follows that
\begin{align*}
|d_t^n u(0)| & = |p_n(f(1), \ldots , d_u^{n-1} f(1))|  \\
          & \le p_n(|f(1)| , \ldots , |d_u^{n-1} f(1)| )\\
          & \le p_n( g(1), \ldots, d_u^{n-1} g(1)) \\
          & = d_t^n v(0),
\end{align*}
demonstrating that $u(t)$ is majorized by $v(t)$ in a neighborhood
of $t=0$.
It follows immediately that
\[
\label{E:uvmajorize}
|u(t)| = \big|\sum_{k=0}^\infty \frac{1}{k!} d_t^{k} u(0) \, t^k \big|
\le \sum_{k=0}^\infty \frac{1}{k!} d_t^{k} v(0) \,  |t|^k \le v(|t|).
\]
The existence of an analytic solution of \eqref{E:MIVODE} with
radius of convergence $|t|<\frac{r}{2 C_\infty}$, given by
\begin{equation}
\label{E:comparison}
v(t) = 1 + r - r\sqrt{1-2C_\infty t / r},
\end{equation}
confirms that $u(t)$ must also be locally analytic about $t=0$.

This argument relies on $C_\infty$, a constant which in practice
is often difficult to ascertain. In our example, it can be
shown that $R=1$, with $r=0.9$, we have
\begin{align*}
C_\infty & = \max_k\{C_0,C_1,C_2, \ldots \} \\
         & = \max_k \big\{
\frac {\mathrm{e}^{2}}{\sin 1}, \frac {9\mathrm{e}^{2}
( 4 \sin 1 - \cos 1 )}{10\sin^2 1}, \frac {81\mathrm{e}^{2}
(2+ 19 \sin^2 1 - 8\cos 1 \sin 1 )}{200 \sin^3 1}, \ldots \big\},
    \end{align*}
and it is not immediately clear where the maximum might occur. An
explicit computation of the $C_k$ terms, plotted in Figure \ref{fig1},
suggests that the maximum occurs near $k=6$,
and one can easily imagine the complexity of $C_k$.
It is also worth noting that the ODE given by  \eqref{E:IVODE}
is simple in
comparison to the systems often used to model
problems of technological importance.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=7cm,height=5cm]{fig1}
\caption{$C_k$ coefficient list} \label{fig1}
\end{center}
\end{figure}

\section{Cauchy solution: the polynomial setting}

We now apply similar techniques to an equivalent polynomial system.
Recall the original problem:
\[
d_t u(t) = \frac{\mathrm{e}^{2 u(t)^2}}{\sin u(t)},
\quad\text{with } u(0) = 1.
\]
Now consider the introduction of the auxiliary variables:
\[
x(t) := \frac{\mathrm{e}^{2 u(t)^2}}{\sin u(t)}, \quad
 y(t) :=\frac{\cos u(t)}{\sin u(t)}.
\]
These auxiliary variables close the system of successive derivatives
that are the foundation of the Cauchy method.  This approach is evident
in the methods of \cite{Griewank,Rall,Jorba2005,Barrio} with an
automatic differentiation flavor, or as suggested by
examples treated in \cite{Perlin1964,Parker1996,Stewart2009}.

We now generate the polynomial system
\begin{gather*}
d_t u  = x   \quad u(0) = 1 \\
d_t x  = \left( 4 x u - x y \right) d_t u  = 4 x^2 u - x^2 y \quad
  x(0) =\frac { \mathrm{e}^{2}}{\sin 1} \\
d_t y  =  -\left(1+y^2 \right) d_t u   = -x-y^2 x  \quad
 y(0) = \cot 1.
\end{gather*}
The first equation is our original ODE; the additional equations
serve a purely computational purpose.

As earlier, we assume the existence of an analytic solution $u$.
We continue by assuming a formal power series for $x$ and $y$,
which can be shown (along with $u$) to be convergent via a
majorant argument.
Now,
\[
u(t) = \sum_{k=0}^\infty u_k t^k , \quad
x(t) = \sum_{k=0}^\infty x_k t^k, \quad
y(t) = \sum_{k=0}^\infty y_k t^k.
\]

The constant on which the previous argument relies is $C_\infty$,
which is difficult in general to construct.
The constants related to the polynomial
argument are easy to construct.
In this new setting, consider the companion problem
\begin{equation}
\label{E:Pmajorize}
d_t z = \mathcal{C} z^m \quad z(0) = c.
\end{equation}
The analytic solution to \eqref{E:Pmajorize} is given explicitly by
\begin{equation}\label{E:MajBound}
z (t) = \left(\mathcal{C}t-\mathcal{C}tm+{c}^{1-m} \right)
^{- (m-1) ^{-1}}.
\end{equation}
In general, \eqref{E:MajBound} will bound solutions to all autonomous
polynomial systems of degree $m$, with suitable choice of $\mathcal{C}$
and initial condition $c$. If $\mathcal{C} = 5, m = 3$ and
 $c =  { \mathrm{e}^{2}}/{\sin 1}$, we
claim that $z(t)$ majorizes $u(t), x(t)$ and $y(t)$.
These parameters arise naturally
when considering the majorization; $\mathcal{C}$ from the largest
row sum of the absolute value of coefficients in the system, $m$
from the largest degree of the polynomial system, and $c$ from
the largest of the absolute value of the initial conditions and 1.
As a brief exercise, we demonstrate this by applying an inductive
argument to verify that the coefficients of the power series
representation of $z(t) = \sum_{k=0}^\infty z_k t^k$ bound those
of $x(t)$.
Clearly $z_0 = c \ge \{|u_0|,|x_0|,|y_0| \}$, since $
c \ge \{ 1, { \mathrm{e}^{2}}/{\sin 1}, \cot 1 \}$
the absolute value of the initial conditions.
Obviously,
\[
z_1 = 5 z_0^3 \ge |4 x_0^2 u_0 - x_0^2 y_0| = |x_1|.
\]
Assuming $z_k \ge \{|u_k|,|x_k|,|y_k| \}$ for $k = 0, \ldots, n$,
it follows that
\begin{align}
z_{n+1} & = \frac{1}{n+1} \cdot 5 \sum_{k=0}^n
\Big( \sum_{i=0}^k  z_i z_{k-i}\Big) z_{n-k} \notag  \\
\label{E:xrecursion}
& \ge \frac{1}{n+1} \cdot \Big| 4 \sum_{k=0}^n
\Big( \sum_{i=0}^k  x_i x_{k-i}\Big) u_{n-k}
- \sum_{k=0}^n \Big( \sum_{i=0}^k  x_i x_{k-i}\Big) y_{n-k} \Big| \\
& = | x_{n+1}| \notag
\end{align}
where a Cauchy product of two series has been applied twice.
An important (and obvious) observation used in \eqref{E:xrecursion}
is that
\[
 x_{n+1}
=\frac{1}{n+1} \cdot \Big[ 4 \sum_{k=0}^n
\Big( \sum_{i=0}^k  x_i x_{k-i}\Big) u_{n-k}
 - \sum_{k=0}^n \Big( \sum_{i=0}^k  x_i x_{k-i}\Big) y_{n-k} \Big],
\]
which can easily be implemented to construct the coefficient
$x_{n+1}$ using
only coefficients of order $n$ or less.
The software tools \texttt{ATOMFT}, \texttt{Taylor}, and most
recently \texttt{TIDES} are
three such packages that exploit this recursive feature
\cite{ATOM,Jorba2005,Abad2009}, although only the last appears
to still be supported.
The polynomial used to construct coefficients in the classic setting,
$p_n$, has now been replaced by an algebraic expression whose complexity
is only $\mathcal{O}(n^3)$. (In fact, augmenting the system allows
reduction to $\mathcal{O}(n^2)$ \cite{Carothers2005}.)
Since $z(t)$ converges on some open interval containing $t=0$
and majorizes $x(t)$ for $|t|<1$, $x(t)$ must also converge on the
intersection of these intervals.
The demonstration is now complete; an explicit verification that
$x(t)$ converges via a term-by-term comparison with the convergent
series representation of $z(t)$.
It is easy to see that a similar argument may be used for $u(t)$
and $y(t)$.

In addition to a simple coefficient recursion and explicit majorization, the
polynomial comparison solution gives rise to an easily computable
local \textit{a priori} error bound. To accomplish this, the
comparison solution $z(t)$ is bounded by $w(t)$, a function with
a geometric series  representation.
We begin with the recurrence relation (see \cite{Warne2006} for a
detailed development) for the coefficients of $z,$
\begin{equation}
\label{E:zrecurrence}
z_{n+1} = \frac{(1+(m-1)n) c^{m-1} \mathcal{C} }{n+1} z_n
\quad z_0 = c, \quad\text{for } n\ge 1.
\end{equation}
For $m \ge 2$,
\begin{equation}\label{E:zbound}
\frac{(1+(m-1)n) c^{m-1} \mathcal{C} }{n+1}
\le (m-1)c^{m-1} \mathcal{C} := \mathcal{C}_\infty.
\end{equation}
Combining \eqref{E:zrecurrence} and \eqref{E:zbound} yields
$z_{n+1} \le \mathcal{C}_\infty z_n$. If
\begin{equation}
\label{E:wrecurrence}
w_{n+1} = \mathcal{C}_\infty w_n, \quad \text{with } w_0 = c,
\end{equation}
then the coefficients of $w$ majorize those of $z$ (and therefore $u$),
and $w(t)$ majorizes $z(t)$ (and $u(t)$). The recurrence relation
\eqref{E:wrecurrence} leads directly to the geometric series,
\[
w(t) = \frac{c}{1-\mathcal{C}_\infty t}
= c \sum_{k=0} (\mathcal{C}_\infty t )^k, \quad\text{when }
|t|<\frac{1}{\mathcal{C}_\infty}.
\]
The function $w$ may be interpreted as a solution to the IVODE
\begin{equation}
\label{E:w}
d_t w(t) = \mathcal{C}_\infty w, \quad w(0) = c
\end{equation}
where $\mathcal{C}_\infty$ bounds the coefficient growth of terms of
$z$, playing much the same role as $C_\infty$.
Here, however, $\mathcal{C}_\infty$
is trivial to compute from \eqref{E:zbound}.

Finally, a simple bound on the remainder term $\mathcal{R}_n$,
given by
\begin{equation} \label{E:errorb}
\mathcal{R}_n(t) := \big| u(t) - \sum_{k=0}^n u_k t^k \big|
 \le c \sum_{k=n+1}^\infty \mathcal{C}_\infty |t|^k
 \le c |\mathcal{C}_\infty t|^{n+1}
\frac{1}{1-|\mathcal{C}_\infty t|},
\end{equation}
provides a concise and computable error bound. For \eqref{E:IVODE},
\[
\mathcal{R}_n(t) \le   \frac {10 \mathrm{e}^{4}}{\sin^2 1}|t|^{n+1}
\frac{1}{1-\frac {10 \mathrm{e}^{4}}{\sin^2 1}|t|}.
\]
For a detailed discussion, and an example for which this bound is tight, see \cite{Warne2006}.
See \cite{Nedialkov} for a detailed discussion of Interval Analysis,
an alternative approach.

In practice, this error bound may be used in a variety of ways.
Here, it suggests a small interval of convergence, with
$t<1/\mathcal{C}_\infty={\sin^2 1}/{10 \mathrm{e}^{4}} \approx 1.3E-3$,
and points to possible singularities in the solution. It can be shown
that the solution of \eqref{E:IVODE} becomes singular quickly
beyond $t \approx 2.6E-2$. It may used to construct a robust
marching method, and provide solutions with known error in regions
of particular interest.

\subsection*{Conclusion}

We have demonstrated that recasting the original ODE as a polynomial
system has several surprising benefits. While the simple differential equation studied here is not tied to a particular modeling scenario, it demonstrates how the conversion makes typically abstract analysis very concrete.
The techniques of Cauchy-Kowalevski, when applied to a polynomial system, lead one naturally to a method in which;
\textit{(i)} coefficients are easily computed in a recursive manner;
i.e., $u_{n+1},x_{n+1}, \text{ and } y_{n+1}$ only depend on products and sums of $\{u_k,x_k,y_k\}_{k=1..n},$
\textit{(ii)} the majorization is explicit, and
\textit{(iii)} there is a clear \textit{a priori} error bound.
The majorization and error bound depend only on immediately observable
quantities of the recast system; coefficient sums, initial conditions,
and degree.

\begin{thebibliography}{10}

\bibitem{Abad2009} A. Abad, R Barrio, F. Blesa, M. Rodriguez;
\emph{a Taylor integrator for differential equations},
 \texttt{http://gme.unizar.es/software/tides},2009.

\bibitem{Barrio} R. Barrio;
 \emph{Performance of the {T}aylor series methods for
{ODEs/DAEs}},
  Appl. Math. and Comput. \textbf{163} (2005), 525--545.

\bibitem{Cauchy} A.-L. Cauchy;
 \emph{M\`{e}moire sur l'emploi du calcul des limites dans
  l'int\'{e}gration des \'{e}quations aux d\`{e}riv\`{e}es partielles}, Comptes
  rendus \textbf{XV} (1842), 44--59.

\bibitem{ATOM} Y. Chang, G. Corliss;
 \emph{\texttt{ATOMFT}: solving {ODEs} and {DAEs} using
  Taylor series}, Computers and Mathematics with Applications \textbf{28}
  (1994), 209--233.

\bibitem{Driver} Bruce K. Driver;
\emph{Cauchy-Kovalevskaya Theorem},
http://www.math.ucsd.edu/~bdriver/231-02-03/Lecture\_Notes/pde4.pdf,
2003.

\bibitem{Evans} L. C. Evans;
 \emph{Partial differential equations}, AMS, Providence, R.I., 1998.

\bibitem{Fehlberg1964} E. Fehlberg;
 \emph{Numerical integration of differential equations by power
  series expansions, illustrated by physical examples}, Tech. Report
  NASA-TN-D-2356, NASA, 1964.

\bibitem{Folland} G. B. Folland;
 \emph{Introduction of partial differential equations}, 2nd ed.,
  Princeton Univ. Press, 1995.

\bibitem{Griewank} A. Griewank, G. F. Corliss;
\emph{Automatic differentiation of algorithms:
 Theory, implementation and application}, SIAM, Philadelphia, PA, 1991.

\bibitem{Jorba2005} A. Jorba, M. Zou;
\emph{A software package for the numerical integration of
  {ODEs} by means of high-order {T}aylor methods},
Exper. Math. \textbf{14} (2005), 99--117.

\bibitem{Kowalevski} S. Kowalevsky;
\emph{Zur theorie der partiellen differentialgleichungen}, J.
  fur reine und angewandte Mathematik \textbf{80} (1875), 1--32.

\bibitem{Nedialkov} N. S. Nedialkov, K. R. Jackson, G. F. Corliss;
\emph{Validated solutions of  initial value problems for ordinary
differential equations}, Appl. Math. and   Comp. \textbf{105} (1999),
21--68.

\bibitem{Parker1996} G. E. Parker, J. S. Sochacki;
\emph{Implementing the {P}icard iteration},
 Neural, Parallel, and Scientific Computation \textbf{4} (1996), 97--112.

\bibitem{Perlin1964} I. E. Perlin, C. P. Reed;
\emph{The application of a numerical integation
  procedure developed by Erwin Fehlberg to the restricted
problem of three  bodies}, Tech. Report NAS8-11129, NASA, 1964.

\bibitem{Pruett2003} C. D. Pruett, J. W. Rudmin,  J. M. Lacy;
\emph{An adaptive {N}-body algorithm of optimal order},
J. Comput. Phys. \textbf{187} (2003), no. 1, 298--317.

\bibitem{Rall} L. B. Rall;
\emph{Automatic differentiation: Techniques and applications},
Lecture Notes in Computer Science, vol. 120, Springer, Berlin, 1981.

\bibitem{Stewart2009} R. D. Stewart, W. Bair;
\emph{Spiking neural network simulation: numerical
  integration with the {P}arker-{S}ochacki method}, J. Comput. Neurosci.
  \textbf{27} (2009), 115--133.

\bibitem{Carothers2005}
D. C. Carothers, G. E. Parker, J. S. Sochacki, P. G. Warne;
\emph{Some properties of solutions to polynomial
  systems of differential equations}, Electron. J. of Diff. Eqns.
\textbf{2005}, no. 40  (2005), 1--18.

\bibitem{Warne2006} P. G. Warne, D. A. Palignone Warne, J. S. Sochacki,
G. E. Parker, D. C. Carothers;
 \emph{Explicit a-priori error bounds and adaptive error control
for approximation of nonlinear initial value differential systems},
Comput. Math. Appl. \textbf{52} (2006), no. 12, 1695--1710.

\bibitem{Paullet} J. E. Paullet, D. A. Polignone Warne;
\emph{Existence and A Priori Bounds
for the Finite Torsion Solution for a Class of Blatz-Ko
Materials}, Math. Mech. Solids \textbf{1} (1996), no. 3, 315--326.

\end{thebibliography}

\end{document}