\documentclass[reqno]{amsart} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2005(2005), No. 40, pp. 1--17.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2005 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2005/40\hfil Some properties of solutions] {Some properties of solutions to polynomial systems of differential equations} \author[D. C. Carothers , G. E. Parker, J. S. Sochacki, P. G. Warne\hfil EJDE-2005/40\hfilneg] {David C. Carothers, G. Edgar Parker, \\ James S. Sochacki, Paul G. Warne} % in alphabetical order \address{Department of Mathematics and Statistics\\ James Madison University\\ Harrisonburg, VA 22807, USA} \email{carothdc@jmu.edu} \email{parkerge@jmu.edu} \email{jim@math.jmu.edu} \email{warnepg@math.jmu.edu} \date{} \thanks{Submitted August 26, 2004. Published April 5, 2005.} \subjclass[2000]{34A05, 34A25, 34A34, 65L05} \keywords{Analytic functions; inverse functions; Maclaurin polynomials; \hfill\break\indent Pade expansions; Gr\"{o}bner bases} \begin{abstract} In \cite{ps1} and \cite{ps2}, Parker and Sochacki considered iterative methods for computing the power series solution to ${\bf y' = G \circ y}$ where ${\bf G}$ is a polynomial from $\mathbb{R}^n$ to $\mathbb{R}^n$, including truncations of Picard iteration. The authors demonstrated that many ODE's may be transformed into computationally feasible polynomial problems of this type, and the methods generalize to a broad class of initial value PDE's. In this paper we show that the subset of the real analytic functions $\mathcal{A}$ consisting of functions that are components of the solution to polynomial differential equations is a proper subset of $\mathcal{A}$ and that it shares the field and near-field structure of $\mathcal{A}$, thus making it a proper sub-algebra. Consequences of the algebraic structure are investigated. Using these results we show that the Maclaurin or Taylor series can be generated algebraically for a large class of functions. This finding can be used to generate efficient numerical methods of arbitrary order (accuracy) for initial value ordinary differential equations. Examples to indicate these techniques are presented. Future advances in numerical solutions to initial value ordinary differential equations are indicated. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem} \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \section{Introduction} ${\bf G}:\mathbb{R}^n \to \mathbb{R}^n$ is an {\em $\mathbb{R}^n$-polynomial} if each component of ${\bf G}$ is a polynomial functional on $\mathbb{R}^n$. A {\em polynomial system} is an autonomous initial value system of differential equations of the form $$ {\bf y}' = {\bf G} \circ {\bf y}, \quad {\bf y}(a)={\bf b} $$ for some $\mathbb{R}^n$-polynomial ${\bf G}$, $a \in \mathbb{R}$, and ${\bf b} \in \mathbb{R}^n$. Let $\mathcal{A}$ be the set of real-analytic functions. An element $f$ of $\mathcal{A}$ is said to be {\em projectively polynomial} if $f$ is a component of the solution to a polynomial system. We denote the collection of all projectively polynomial functions by $\mathcal{P}$. $\mathcal{P}^a$ will denote those elements of $\mathcal{P}$ with initial condition at $a \in \mathbb{R}$. In \cite{ps1} and \cite{ps2}, Parker and Sochacki showed that many differential equations may be re-written as polynomial systems, and also gave efficient algorithms for the numerical solution of these systems. In this paper, we further explore the properties of solutions to such systems, and we begin by classifying elements of $\mathcal{P}$ according to the dimension of the polynomial system and degree of the underlying polynomial for which they are components of the solution. We show in Theorem 1 that every element of $\mathcal{P}$ is the solution to a system of at most second-degree. In Theorem 2, we show that there are real-analytic functions that are not elements of $\mathcal{P}$ and establish that $\mathcal{P}$ shares the field and near-field structure of the real-analytic functions, making it a proper subalgebra in both senses. In Theorem 3, we use Gr\"{o}bner basis techniques to show that membership of a particular function in $\mathcal{P}$ is equivalent to that function being a solution to a single higher order polynomial differential equation in a single variable. Example 5 will use this result to show how a wide range of systems of differential equations may be ``de-coupled" to produce equivalent equations in a single variable. Note that elementary facts about differential equations show that $\mathcal{P}$ contains the functions sine, cosine, exponential, and polynomials. In \cite{ps1} and \cite{ps2}, it was shown that $\mathcal{P}$ is closed under addition, multiplication, and function composition. Theorem 4 below shows that local inverses of elements of $\mathcal{P}$ are also in $\mathcal{P}$. Function $r$ defined by $r(t) = {{1}\over{t}}$ is the solution to $$ r' = -r^2; \quad r(1) = 1 $$ and thus $r \in \mathcal{P}$. If $h(t) = t^p$ for $p \in \mathbb{R}$, we note that $h$ and $r$ are the solutions to the system \begin{gather*} h' = phr \\ r' = -r^2\\ h(1) = 1, \quad r(1) = 1 \end{gather*} so that $h \in \mathcal{P}$. Hence, $\mathcal{P}$ includes all algebraic combinations and local inverses of the elementary functions. As an illustration of how large is the class of differential equations that may be transformed into polynomial form, Theorem 5 will demonstrate that if $f \in \mathcal{P}$ and $y$ is a solution to the differential equation $y' = f \circ y$, then $y \in \mathcal{P}$. Finally, in Theorem 6 we demonstrate cosets for $\mathcal{P}$ by $*$ and $\circ$ can intersect only in elements of $\mathcal{P}$ except when represented by the respective identity elements. Applications that are consequences of the algebraic structure are suggested. We use the development discussed in the last paragraph to show that we can use Cauchy's theory for the analytic solution to an initial value ordinary differential equation with analytic generator and Cauchy products to improve on the methods presented in \cite{ps1} and \cite{ps2}. We give examples to show how the theory presented can be used in a numerical or symbolic environment. As a matter of fact, one can build a computing environment that does strictly polynomial arithmetic to any order of accuracy desired. The need to have Maclaurin or Taylor polynomials stored in a computer library is only necessary for the initial conditions since we can generate the Maclaurin and Taylor polynomials for the solution to an initial value ordinary differential equation with polynomial right hand side to any degree desired with one algorithm. \section{$\mathcal{P}$ and Polynomial Degree} For $f \in \mathcal{P}$, if $f$ is a component of the solution to ${\bf y}' = {\bf G} \circ {\bf y}$, ${\bf y}(a)={\bf b}$, then ${\bf G}$ will be called a {\em generator} for $f$. $\mathcal{P}_{i,k} \subseteq \mathcal{P}$ will denote those elements having a polynomial generator whose domain has dimension that is less than or equal to $i$ and whose (polynomial) degree is less than or equal to $k$. $\mathcal{P}_k$ will denote $\cup_{i\in \mathbb{N}} \mathcal{P}_{i,k}$. For $f \in \mathcal{P}$, a particular system of differential equations for which $f$ is a component of the solution will be called a {\em projection} for $f$. In \cite{ps1} and \cite{ps2} examples were given that indicate that lower degree, higher dimension projections may be more numerically robust than related lower dimension, higher degree projections. In \cite{ps1}, a first degree projection was shown for each polynomial, and $\mathcal{P}_1$ is thus dense in $\mathcal{P}$. The following theorem shows that every element of $\mathcal{P}$ has a second-degree projection. \begin{theorem} $\mathcal{P} = \mathcal{P}_2$. \end{theorem} \begin{proof} Suppose $f \in \mathcal{P}$, and that $f$ is the first component of a solution to $$ {\bf y}'={\bf G} \circ {\bf y} = {\bf G}(y_1,\dots ,y_n), ~{\bf y}(a)={\bf b}. $$ If $\{G_1,\dots ,G_n\}$ is the set of components of ${\bf G}$, let $N_j$ denote the largest power of $y_j$ that appears in any of the polynomials $\{G_1,\dots ,G_n\}$. We make the substitution $$ v_{i_1,\dots ,i_n} = y_1^{i_1}y_2^{i_2}\dots y_n^{i_n} $$ with each $0 \leq i_j \leq N_j$ and at least one of $i_j>0$. Note that under this substitution $y_1 = v_{1,0,\dots ,0}, y_2 = v_{0,1,0,\dots ,0}$, and so on. Substituting into each of the original equations for $y_j'$ yields a first degree equation in the variables $\{v_{i_0,\dots ,i_n}\}$, in which the left side is $v'_{i_0,i_1,\dots ,i_n}$ with only one $i_j = 1$ and the remaining $i_j = 0$. This set of equations may be augmented with additional equations of at most second degree to form a new polynomial system by adding an additional equation for each of the remaining $v_{i_1,\dots ,i_n}$ variables that does not correspond to one of the original $y_j$. For these variables, we have formally $$ v_{i_1,\dots ,i_n}'= \sum_{k=0}^n i_k v_{i_1,\dots ,i_k-1,\dots ,i_n} y_k'. $$ Since each $y_k'$ may be replaced under the substitution by a linear combination of the variables $\{v_{i_0,\dots ,i_n}\}$, each equation in the new system has degree at most two. $f$ is the first component (corresponding to $v_{1,0,\dots,0}$) of a solution to this system. \end{proof} \section{Characterization of $\mathcal{P}$} Recall that $\mathcal{P}^a$ denotes those elements of $\mathcal{P}$ with initial condition at $a \in \mathbb{R}$. In \cite{ps1}, $\mathcal{P}^a$ is shown to be closed under addition and multiplication, and that $f \circ g \in \mathcal{P}^a$ for $g \in \mathcal{P}^a$ and $f \in \mathcal{P}^{g(a)}$. For any fixed real number $\alpha$, suppose that $h$ is defined by $$ h(t) = {{1}\over{t-\alpha}}, \quad t>\alpha. $$ Then $h \in \mathcal{P}^a$ for any $a > \alpha$. This guarantees the local invertibility of any non-zero element of $\mathcal{P}$ under multiplication. Thus $\mathcal{P}$ shares with $\mathcal{A}$ both its local field properties under addition and multiplication and its local near-ring properties under addition and function composition. First we show that $\mathcal{P}$ is a proper subset of $\mathcal{A}$. \begin{theorem} $\mathcal{A} \neq \mathcal{P}$. \end{theorem} \begin{proof} We will find an analytic function $u_1 \notin \mathcal{P}$ by specifying values for $u_1^{(j)}(0)$, $|u_1^{(j)}(0)| \leq 1$ for the series expansion of $u_1$ in such a way that $u_1 \notin \mathcal{P}_{m,2}$ for any $m$. It is possible to choose values for $u_1^{(j)}(0)$, $0 \leq j \leq 4$, $|u_1^{(j)}(0)| \leq 1$, so that any function whose 4-th degree Maclaurin polynomial is $$ \sum_{j\leq 4}{{u_1^{(j)}(0)}\over{j!}} t^j $$ is not in $\mathcal{P}_{1,2}$. To show this, note that $u_1 \in \mathcal{P}_{1,2}$ implies that there are numbers $a_0, a_1, a_2$ so that $$ u_1' = a_0 +a_1u_1 + a_2u_1^2 $$ Differentiate three times: \begin{gather*} u_1'' = a_1u_1' + 2a_2u_1u_1' \\ u_1''' = a_1u_1'' + 2a_2(u_1u_1''+u_1^{\prime 2})\\ u_1^{(4)} = a_1u_1''' + 2a_2(u_1u_1'''+3u_1'u_1'') \end{gather*} Any choice of values for $u_1^{(j)}(0)$, $0 \leq j \leq 4$ yields a system of 4 (linear) equations in the unknowns $a_0, a_1, a_2$. It is of course possible to make this choice so that the system of equations has no solution (e.g. $u_1^{(j)}(0) = 1$ for $j= 0, 1, 2, 3$ and $u_1^{(4)}(0)= {{1}\over{2}}$). For any such function $u_1$, there is no appropriate differential equation which $u_1$ satisfies, and we conclude $u_1 \notin \mathcal{P}_{1,2}$. Now assume that for a particular value of $K$, values for $u_1^{(j)}(0)$, $0 \leq j \leq n_K$, $|u_1^{(j)}(0)| \leq 1$, have been determined so that any function whose $n_K$ degree Maclaurin polynomial is $$ \sum_{r\leq n_K}{{u_1^{(r)}(0)}\over{r!}} t^r $$ is not in $\mathcal{P}_{K,2}$. In order that $u_1 \in \mathcal{P}_{K+1,2}$, it would be necessary that there exist real analytic functions $u_2,\dots ,u_{K+1}$ and numbers $a_{ijk}$, $b_{ij}$, $c_i$ with $0 \leq i,j, k \leq K+1$, so that $u_i$, $1 \leq i \leq K+1$ satisfy the system $$ u_i'=\sum_{j=1}^{K+1}\sum_{k=1}^{K+1}a_{ijk}u_ju_k+\sum_{j=1}^{K+1}b_{ij}u_j+c_i. $$ Evaluate $u_i'$ at $0$ to obtain $$ u_i'(0)=Q_{1,i} $$ where each $Q_{1,i}$ is a polynomial in the unknowns $\{a_{ijk}, b_{ij}, c_i, u_2(0),\dots ,u_{K+1}(0)\}$ (let $M$ represent the number of these unknowns). Differentiate each of the $u_i'$, evaluate at $0$, and replace each $u_i'(0)$, $2 \leq i \leq K+1$ with the corresponding $Q_{1,i}$ to obtain $$ u_i''(0) = Q_{2,i} $$ where each $Q_{2,i}$ is a polynomial in the same set of $M$ unknowns. Proceeding in this way, we obtain a set of polynomial equations $$ u_i^{(r)}(0)=Q_{r,i}\,,\quad 1\leq r \leq n_K $$ for $1 \leq i \leq K+1$ (note that the left side of the equations with $i=1$ are already defined). We will consider solutions in the $M$ unknowns to the equations whose left side is $u_1^{(j)}(0)$. These equations determine a (possibly empty) solution set $S \subseteq \mathbb{C}^M$. Let $dim(S)$ denote the dimension (as an algebraic variety) of $S$ in $\mathbb{C}^M$. Note that $dim(S)< M$. Differentiate $u_1^{(n_K)}$. Evaluate at $0$ and replace each $u_i^{(n_K)}(0), \quad 2 \leq i \leq K+1$ with $Q_{n_K,i}$ as before to obtain a polynomial equation in the $M$ unknowns $$ u_1^{(n_K+1)}(0) = Q_{n_K+1,1} $$ with $u_1^{(n_K+1)}(0)$ not yet determined. For any choice of $u_1^{(n_K+1)}(0) $, the system of equations $$ u_1^{(r)}(0)=Q_{r,1}\,,\quad 1\leq j\leq n_K+1 $$ determines a solution set $S^* \subseteq S$. $u_1^{(n_K+1)}(0)$ may be chosen so that $dim(S^*) n_K+1 $$ where each new polynomial reduces the dimension of the solution set, to produce choices for $u_1^{(n_K+1)}(0),\dots ,u_1^{(n_{(K+1)})}(0) $ so that the resulting solution set in $\mathbb{C}^M$ is empty. Thus, for any function whose $n_{(K+1)}$ degree Maclaurin polynomial is $$ \sum_{r\leq n_{(K+1)}}{{u_1^{(r)}(0)}\over{r!}} t^r $$ it is impossible that $u_1 \in \mathcal{P}_{K+1,2}$. We may proceed as above to inductively construct $u_1$ so that $u_1 \notin \mathcal{P}_{n,2}$ for any $n$, noting that $|u_1^{(r)}(0)| \leq 1$ for all $r$ guarantees the convergence of the series. \end{proof} Next, we provide an alternative characterization of the elements of $\mathcal{P}$. \footnote{The authors thank Professor Carter Lyons for his contributions to this theorem.} One consequence of the following theorem is the possibility of ``de-coupling'' any system of differential equations that may be re-cast as a polynomial system, and writing a single differential equation for each of the original variables. (See Example 5.) \begin{theorem} \label{thm3} $u \in \mathcal{P}$ if and only if for some $n$ there is a polynomial $Q$ in $n+1$ variables so that $Q(u, u', \dots, u^{(n)})=0$. \end{theorem} \begin{proof} ($\Rightarrow$) If $u$ is projectively polynomial, then there are functions $v_1, \dots, v_{n-1}$ and polynomials $P_0, \dots P_{n-1}$ so that \begin{gather*} u' = P_0(u, v_1, \dots, v_{n-1}),\\ v_j' = P_j(u, v_1, \dots, v_{n-1}) \end{gather*} for $j=1,\dots,n-1$. Formally differentiate both sides of the first equation to find $u''$ and substitute for $u',\dots, v_{n-1}'$ in the right side to write $u''$ as a polynomial in terms of $u, \dots,v_{n-1}$. Repeat this procedure several times on the result to find $u''',\dots,u^{(n)}$, obtaining: $$ u^{(j)}=P_{j}^*(u, v_1,\dots, v_{n-1}) $$ for $j=2,\dots,n$. For convenience we will rename $P_0=P_1^*$ and let \\ $P_0^*(u, v_1,\dots, v_{n-1}) =u$ so that $u^{(j)} = P_j^*(u, v_1,\dots, v_{n-1})$ for $j=0,\dots,n$. For the moment we will treat the ``$u$'' and ``$v$'' functions as strictly algebraic objects, and the above as polynomial equations in $2n+1$ variables: $$ \{v_1,\dots,v_{n-1},u,u^{(0)},u',\dots,u^{(n)}\} $$ We distinguish between $u$ and $u^{(0)}$ (as algebraic objects), and will assume this order of the variables in what follows. We seek a polynomial differential equation in the last $n+1$ variables $$ \{ u^{(0)}, u', \dots, u^{(n)}\} $$ in which the equation is satisfied by all points given by a solution to the original system of differential equations. To this end, we wish to algebraically manipulate the above set of equations to eliminate the variables $\{v_1,\dots,v_{n-1},u\}$. Let $I$ denote the ideal in the set of all polynomials in $\{v_1,\dots,v_{n-1},u,u^{(0)},u',\dots,u^{(n)}\}$ generated by the collection $$ \{u^{(j)} - P_{j}^*, j=0,\dots,n\}, $$ and let $I_{n}$ be those elements of $I$ that involve only the last $n+1$ variables \\ $\{u^{(0)}, u', \dots, u^{(n)}\}$. $I_{n}$ is referred to as the the $n$th {\em elimination ideal}. The equations $\{u^{(j)} = P_j^*(u, v_1,\dots, v_{n-1}), j=0, \dots, n \}$ are a {\em parameterization} for an {\em algebraic variety} $V$ in $\mathbb{C}^{n+1}$. Here the term {\em variety} denotes a collection of points in $\mathbb{C}^{n+1}$ determined by a set of polynomial equations in $\{u^{(j)}\}$. The points in $\mathbb{C}^{n+1}$ with coordinates $u^{(j)}$ determined by allowing the {\em parameters} $\{v_1,\dots, v_{n-1},u\}$ to vary in $\mathbb{C}^{n}$ generate (in an alternate way) a variety for some set of polynomial equations in $\{u^{(j)}\}$, technically the smallest variety containing the points given by the parameterization, in case these points do not fill up the smallest variety containing them. In the language of algebraic geometry, we seek an {\em implicitization} that will provide the polynomial equations in $\{u^{(j)}\}$ alone defining $V$, with the parameters eliminated. Applying Gr\"{o}bner basis elimination techniques with a lexicographic ordering determined by our order of the variables above, we conclude that $I_{n}$ contains non-zero elements. Specifically, by Chapter 3, Section 3, Theorem 1 of \cite{clo}, if $\bf{G}$ is a Gr\"{o}bner basis for $I$, then $I_{n}\cap\bf{G}$ is a Gr\"{o}bner basis for $I_{n}$. The comments in \cite{clo} following this theorem describe the specifics of the algorithm for eliminating the variables $\{v_1,\dots,v_{n-1},u\}$. A basis element $Q \in I_{n}\cap\bf{G}$ gives us the appropriate differential equation $Q(u^{(0)}, u', \dots, u^{(n)})=0$ (and at this point we may for convenience identify $u^{(0)}$ with $u$). In this way, the resulting Gr\"{o}bner basis provides an equation in $\{u^{(0)}, u', \dots, u^{(n)}\}$ that includes in its solution set any solution to the original (parameterized) equations $\{u^{(j)} = P_j^*(u, v_1,\dots, v_{n-1})\}$ \smallskip To prove ($\Leftarrow$), suppose that $u$ is a solution to $$ Q(u, u', \dots, u^{(n)})=0 $$ for some polynomial $Q$. We introduce new variables by letting $$ v_j=u^{(j)} $$ for $j=0,\dots,n$, noting that $u=v_0$. The above polynomial equation becomes $$Q(v_0, v_1, \dots, v_n)=0.$$ Note that for each $i$, ${{\partial Q}\over{\partial v_i}}$ is a polynomial in the collection $\{v_j\}$, as are the second partials. We introduce one additional variable $w$ corresponding to $1/{{\partial Q}\over{\partial v_n}}$. The desired polynomial system is then \begin{gather*} v_j' = v_{j+1},\quad j=0,\dots,n-1 \\ v_n' = -w\sum_{j=0}^{n-1}{{\partial Q}\over{\partial v_j}}v_{j+1} \\ w' = -w^2(\sum_{j=0}^{n-1}{{\partial^2 Q}\over{\partial v_j\partial v_n}}v_{j+1}+{{\partial^2 Q}\over{\partial v_n^2}}(-w\sum_{j=0}^{n-1}{{\partial Q}\over{\partial v_j}}v_{j+1}) ). \end{gather*} \end{proof} The following theorem shows that $\mathcal{P}$ is a local sub-near-field of $\mathcal{A}$ by showing that local inverses under composition of elements of $\mathcal{P}$ are elements of $\mathcal{P}$. \begin{theorem} \label{thm4} Suppose $f \in \mathcal{P}^a$, $f'(a) \neq 0$ and $f^{-1}$ is a local inverse for $f$ in a neighborhood of $a$. Then $f^{-1} \in \mathcal{P}^{f(a)}$. \end{theorem} \begin{proof} Let $f \in \mathcal{P}^a$ with local inverse $f^{-1}$ in a neighborhood of $a$. There is an open interval $U$ over which $f$ has a non-zero derivative. Since $f \in \mathcal{P}^a$, we have $f \in \mathcal{P}^a_{n,2}$ for some $n$ by Theorem 1. For convenience, we let $f_1=f$, and note that there are functions $f_i, ~~2 \leq i \leq n$ and numbers $a_{ijk}$, $b_{ij}$, $c_i$ with $1 \leq i,j, k \leq n$, so that $$ f_i'=\sum_{j=1}^{n}\sum_{k=1}^{n}a_{ijk}f_jf_k+\sum_{j=1}^{n}b_{ij}f_j+c_i. $$ Note that if we differentiate the above and replace every resulting occurrence of $f_j', f_k'$ in the right side with a corresponding polynomial in $f_j, f_k$, we may note that $f_i''$ is a polynomial in the original functions. If $t$ is an element of the range of $f$ restricted to $U$, then $$ f'\circ f^{-1}(t)(f^{-1})'(t) = 1 $$ or equivalently $$ (f^{-1})'(t) = {{1}\over{f'\circ f^{-1}(t)}}. $$ Also note that $$ ({{1}\over{f'\circ f^{-1}}})' = -({{1}\over{f'\circ f^{-1}}})^2 ((f''\circ f^{-1})(f^{-1})')=-({{1}\over{f'\circ f^{-1}}})^3(f''\circ f^{-1}). $$ In the last expression, it is possible to replace $(f''\circ f^{-1})$ with a polynomial in the collection $\{f_i \circ f^{-1}\}$. We now consider a system of polynomial equations in $f_1^{-1}$, ${{1}\over{f' \circ f_1^{-1}}}$, and $\{f_i \circ f^{-1}, 1 \leq i \leq n\}$. From the above, we see that $(f^{-1})'$ and ${{1}\over{f' \circ f_1^{-1}}}'$ may be expressed as polynomials in this collection of functions. In addition, $$ (f_i \circ f^{-1})' = (f_i' \circ f^{-1})(f^{-1})' $$ may be expressed as a polynomial in this collection of functions by replacing $f_i'$ with $\sum_{j=1}^{n}\sum_{k=1}^{n}a_{ijk}f_jf_k+\sum_{j=1}^{n}b_{ij}f_j+c_i$ \end{proof} The above theorem guarantees that $\mathcal{P}$, like $\mathcal{A}$, is a local near-field for $+$ and $\circ$. In Section 5, we will see that reasonably tractable projections can occur when considering certain inverse problems. A result of Cauchy (\cite{br}, p. 169, Corollary 2) shows that if $G$ is real-analytic then the solution to $$ y' = G \circ y $$ must also be real-analytic on some open interval containing its initial condition. Note that a computation similar to that in the proof of the preceding theorem shows that if $f$ is real-analytic, $f'$ is not 0 on an open interval and $g$ is the local inverse of $f$ then $$ f'(t) = \frac{1}{g'(f(t))}. $$ That is, the converse of Cauchy's theorem is true. \footnote{The authors thank Professor Peter Lax for his insight on this observation.} >From classical analysis if $f$ is a real-analytic function, $y_0$ is interior to interval of convergence of $f$, and $y$ is a solution to $$ y' = f \circ y; \quad y(0) = y_0 $$ then $y$ is also real analytic. The following theorem shows that an analogous statement is also true with ``real analytic'' replaced by ``projectively polynomial,'' demonstrating something of the range of non-polynomial differential equations that may be transformed into polynomial systems. \begin{theorem} \label{thm5} Suppose that $f \in \mathcal{P}^{y_0}$. If $y$ is a solution to $$y' = f \circ y; ~~~y(0)= y_0$$ then $y \in \mathcal{P}^0$. \end{theorem} \begin{proof} Since $f \in \mathcal{P}^{y_0}$ we may assume there are $n$ $\mathbb{R}^n$-polynomials $Q_j$ so that that $f$ is the first component (setting $f_1=f$) of the solution to the system $$ f_j' = Q_j(f_1,\dots,f_n); \quad f_j(y_0)=b_j, \quad j= 1, \dots, n. $$ Thus \begin{align*} (f_j\circ y)' &= (f_j' \circ y) y' \\ &= (Q_j(f_1,\dots,f_n)\circ y)(f_1\circ y) \\ &= Q_j(f_1\circ y,\dots,f_n \circ y)(f_1\circ y). \end{align*} The above along with the original equation $y' = f_1 \circ y$ provides a polynomial system for which the components of the solution are $y$ and $f_j \circ y$, $ j = 1, \dots, n$. Hence $y \in \mathcal{P}^0$. \end{proof} \section{The intersecting coset property} The previous section establishes that $\mathcal{P}$ is a proper sub-algebra of $\mathcal{A}$ both as a local field by $+$ and $*$ and as a local near-field by $+$ and $\circ$. In this section we examine consequences of the fact that $\mathcal{P} - \{ 0 \}$ under $*$ and the non-constant elements of $\mathcal{P}$ under $\circ$ are local groups. \begin{theorem} \label{thm6} Suppose that $f \in C^1$ is defined in a neighborhood of $a \in \mathbb{R}$ and that $\mathcal{P}^{f(a)}\circ f \cap f*\mathcal{P}^a \neq \{kf | k \in \mathbb{R}\}$. Then $f \in \mathcal{P}$. \end{theorem} \begin{proof} Since constant functions are elements of $\mathcal{P}$, suppose that $f$ is not a constant function and that $g \in \mathcal{P}^{f(a)}$ and $ h \in \mathcal{P}^a$ satisfy the following equation $$ g \circ f = h *f, $$ with $h$ not a constant function. Since $g \in \mathcal{P}^{f(a)}$, there are functions $g_j,~ 2 \leq j \leq n$ defined in a neighborhood of $f(a)$ (and we let $g_1=g$ for convenience) so that each $g_i'$ may be written as a polynomial expression in the collection $\{g_j\}$. Since $ h \in \mathcal{P}^a$, we may similarly find a collection $h_j, ~2 \leq j \leq m$ defined in a neighborhood of $a$ (letting $h=h_1$) satisfying a similar system of polynomial differential equations. We restrict attention to an interval containing $a$ that is in the common domain for $f$ and the collections $\{g_j \circ f\}$ and $\{h_j\}$. Because $g \circ f = h *f$, we have $$ (g' \circ f)f'= h' f+ hf'. $$ If $g' \circ f -h$ is not identically zero, let $U$ be an open interval on which $g' \circ f -h$ is non-zero. Then for $t \in U$ we have $$ f' (t)= {{h'(t)f(t)}\over{g' \circ f(t)- h(t)}} . $$ Consider a system of polynomial equations involving $f$, ${{1}\over{g' \circ f- h}}$, $\{g_j \circ f, 1 \leq j \leq n$\}, and $\{h_j, 1 \leq j \leq m\}$. The above shows that $$ f' = h_1' f {{1}\over{g' \circ f- h}} $$ and we may express $f'$ as a polynomial in this collection by replacing $h_1'$ with the corresponding polynomial in the $h_j$. Also $$ ({{1}\over{g' \circ f- h}})' = - ({{1}\over{g' \circ f- h}})^2(((g''\circ f)*f')-h'). $$ Since $g'=g_1'$, $g'$ may be written as a polynomial in the collection $\{g_j\}$. Differentiating both sides and replacing each $g_j'$ on the right side with a polynomial in the collection $\{g_j\}$, we see that $g''$ may be written as a polynomial in the collection $\{g_j\}$. $h' =h_1'$ may be written as a polynomial in the collection of $\{h_j\}$. Thus, the above equation shows that $({{1}\over{g' \circ f- h}})'$ may be written as a polynomial in $f$, ${{1}\over{g' \circ f- h}}$, $\{g_j \circ f | 1 \leq j \leq n$\}, and $\{h_j | 1 \leq j \leq m\}$. Finally, since $$ (g_i \circ f)' = (g_i' \circ f)*f' $$ we note again that $g_i'$ may be replaced by a polynomial in the collection $\{g_j\}$, so that $g_i'$ may also be written as a polynomial in $f$, ${{1}\over{g' \circ f- h}}$, $\{g_j \circ f, 1 \leq j \leq n$\}, and $\{h_j, 1 \leq j \leq m\}$. The above shows that $f \in \mathcal{P}$ as long as $(g' \circ f) -h$ is not identically $0$ (and note that $f \in \mathcal{P}^a$ in case the interval $U$ contains $a$). Otherwise, since in this case $g' \circ f =h$ we have $f= (g')^{-1}\circ h$ and thus $f \in \mathcal{P}^a$, as long $g'$ is not constant. But $g'$ constant would imply that $h$ is constant in the original relation $g \circ f = h*f$. \end{proof} \section{Examples} \noindent {\bf Example 1.} Consider the question of finding a solution to the polynomial equation $\sum_{k=0}^na_k t^k=0$, $a_n \neq 0$. Let $f$ be the polynomial function defined by $$ f(t) = \sum_{k=0}^na_k t^k. $$ >From the proof of Theorem 4 we conclude $$ (f^{-1})' = {{1}\over{f'\circ f^{-1}}} $$ and $$ ({{1}\over{f'\circ f^{-1}}})' = -({{1}\over{f'\circ f^{-1}}})^3(f''\circ f^{-1}). $$ Thus, any local inverse for $f$ is a solution to the system \begin{gather*} x'=y \\ y'=-y^3\sum_{k=2}^nk(k-1)x^{k-2} \end{gather*} with the particular local inverse determined by the initial conditions. If $f'(p)\neq0$ and we use the initial conditions $x(f(p))=0$ and $y(f(p)) ={{1}\over{f'(p)}}$, then computation of $x(0)$ solves the polynomial equation. Notice that the application of the methods of \cite{ps1} to this system provides an alternative to the classical power series methods discussed in \cite{ap}. We have a qualitative guarantee associated with this computation; we will find the zero of $f$ in the portion of the domain containing $p$ over which $f$ is invertible, if there is such a point, or else the second component will increase or decrease without limit. One may compare this with Newton's method, where conditions on the geometry of the graph of the function modelling the equation are necessary to guarantee the convergence to a particular zero of the equation and failure of the iteration is no guarantee of non-existence of a zero within the component of the initial guess. \smallskip According to Theorem 1 of this paper, if we have an algorithm to generate the Maclaurin polynomial for the solution to the initial value problem $$ {\bf y}' = {\bf Q} \circ {\bf y}, ~~{\bf y}(0)={\bf a} $$ where $ {\bf Q} $ is linear or quadratic on $\mathbb{R}^n$, and ${\bf a} \in \mathbb{R}^n$. Then we can generate the Maclaurin polynomial for the solution to any initial value ordinary differential equation (IVODE) with a polynomial right hand side. From Theorem 1 of \cite{ps1} we know that the Picard iterates contain the Maclaurin polynomial for the solution to the initial value ordinary differential equation, and the proof in that paper shows that only that part of the iterate which is the Maclaurin polynomial need be carried forward to continue the process. In what follows, we show that we can further modify the Picard iteration process so that the Picard iterates generate the terms of the Maclaurin polynomial algebraically. We illustrate the process in $\mathbb{R}^1$, taking $ Q(y) = A+By+Cy^2 $ and $ y(0) = \alpha. $ Start the Picard process with $$ P_1(t) = \alpha = M_1(t) $$ and then define $$ P_{k+1}(t) = \alpha + \int_{0}^t Q(P_k(s))ds = \alpha + \int_{0}^t (A + B P_k(s) + C P_k(s)^2)ds. $$ This integral equation is equivalent to $$ P_{k+1}'(t) = A + B P_k(t) + C P_k(t)^2 = Q(P_k(t)) \quad \ y(0) = \alpha. $$ In \cite{ps1} we showed that $$ P_k(t) = M_k(t) + \sum_{n=k}^{2^k-1} b_n t^n, $$ where $$ M_{k+1}(t) = \sum^k_{n=0} a_n t^n $$ is the degree $k$ Maclaurin polynomial for the solution $y$ to the initial value ordinary differential equation. One can use the integral equation or differential equation for $P_{k+1}(t)$ to solve for $a_k$ and $b_k$. Since $a_0 = \alpha$ and $a_1 = Q(\alpha)$, this only needs to be done for $k \geq 2$. Using Cauchy products this leads to \begin{align*} P_{k+1}' &= \sum^k_{n=1} n a_n t^{n-1} + \sum_{n=k+1}^{2^{k+1}-1} n b_n t^{n-1} \\ &= \sum^{k-1}_{n=0} (n+1)a_{n+1} t^n + \sum_{n=k+1}^{2^{k+1}-1} n b_n t^{n-1} \\ &= A + B \sum^{2^k-1}_{n=0} b_n t^n + C \sum^{2^k-1}_{n=0} b_n t^n \sum^{2^k-1}_{n=0} b_n t^n \\ &= A + B \sum^{2^k-1}_{n=0} b_n t^n + C \sum^{2^k-1}_{j=0} \sum^{2^k-1}_{n=0} b_j b_n t^{n+j} \\ &= A + B \sum^{2^k-1}_{n=0} b_n t^n + C \sum^{2(2^k-1)}_{n=0} d_n t^n, \end{align*} where $b_n=a_n$ for $n \leq k$ and $d_n = \sum_{j=0}^n b_j b_{n-j}$ for $0 \leq n \leq 2^k-1$ and $d_n = \sum_{j={n-2^k+1}}^{2^k-1} b_j b_{n-j}$ for $2^k \leq n \leq 2(2^k-1).$ By equating like powers it is straightforward to show that $$ k a_k = B a_{k-1} + C \sum^{k-1}_{j=0} a_j a_{k-1-j}. $$ That is, $$ M_{k+1}(t) = M_k(t) + \frac{(B a_{k-1} + C \sum^{k-1}_{j=0} a_j a_{k-1-j} )}{k} t^k. $$ That is, we can obtain the Maclaurin polynomial algebraically. This is easily extended to any polynomial ODE or system of polynomial ODE's. Therefore, for ODE's with polynomial right hand side the Maclaurin coefficients of the solution can be obtained algebraically with the same algorithm. This algorithm we call the Algebraic-Maclaurin algorithm. The Cauchy result \cite[p. 169, Corollary 2]{br} mentioned after the proof of Theorem 4 shows that this algorithm gives the analytic solution to the initial value ordinary differential equation. One can also generate a formula for the $b_n's$ for $n >k$ using the above results. In a future paper, we exploit this to give convergence rates and error estimates for polynomial differential equations. It is also easy to see that one can modify the above algorithm to work for any polynomial system. In the examples below we do this in our numerical environments. We now present three examples using the Algebraic-Maclaurin algorithm. Since this algorithm is totally algebraic, we can generate as many Maclaurin coefficients as wanted at any time step in the algorithm. Therefore, we have a method that can change the order (degree) of accuracy on the 'fly'. We compare numerical results from this method with numerical results from the fourth order Runge-Kutta method. \noindent{\bf Example 2.} Consider the initial value ordinary differential equation (I): \begin{equation} \label{I} u' = -u^r ; \; u(0) = \alpha. \end{equation} If we let $x_1=u$, $x_2=u^r$ and $x_3=x_1^{-1}$, we obtain (II): \begin{equation} \label{II} \begin{gathered} x_1' = -x_2; \quad x_1(0)= \alpha \\ x_2' = -rx_2^2x_3;\quad x_2(0) = \alpha^r \\ x_3' = x_2x_3^2; \quad x_3(0)= \alpha^{-1}. \end{gathered} \end{equation} If we let $x_4=x_2x_3$, we obtain (III): \begin{equation} \label{III} \begin{gathered} x_1'= -x_2; \quad x_1(0)= \alpha \\ x_2'= -rx_2x_4;\quad x_2(0) = \alpha^r \\ x_3' = x_3x_4; \quad x_3(0)= \alpha^{-1} \\ x_4' = (1-r)x_4^2; \quad x_4(0)=\alpha^{r-1}. \end{gathered} \end{equation} Noticing from this that we do not need the $x_3$ equation, we obtain (IV): \begin{equation} \label{IV} \begin{gathered} x_1' = -x_1x_4; \quad x_1(0) = \alpha \\ x_4' = (1-r)x_4^2; \quad x_4(0) = \alpha^{r-1}. \end{gathered} \end{equation} Since the solution for the original initial value ordinary differential equation can be obtained from this last IVODE, $x_1$ is in $\mathcal{P}_{2,2}$. We compare the solutions to \eqref{I}--\eqref{IV} using the fourth order Runge-Kutta scheme and the Algebraic-Maclaurin algorithm to generate the fourth, fifth and sixth order (degree) Maclaurin polynomials for \eqref{II}--\eqref{IV}. Of course, the Maclaurin polynomial for $u$ is the same in \eqref{I}--\eqref{IV}. However, in a numeric environment we update in $t$ using numerical values which leads to different numerical values for each method. In the table below is the error of the results of a simulation with a time step of 0.0625 using fourth order Runge-Kutta on the systems \eqref{I}--\eqref{IV} and the fourth, fifth and sixth order Maclaurin polynomial on \eqref{II}--\eqref{IV} using the algorithm above with $r=3/2$ and $x(0)=1.25$. The computing times were essentially the same. (See Example 4 below where computing times are relevant.) The errors are in comparison with the exact solution given by Maple using 30 digits of accuracy. \begin{center} \begin{tabular}{|c|c|c|c||c|} \hline Time & R-K on I & R-K on II & R-K on III & R-K on IV \\ \hline\hline $0.125$ & 7.03864E-08 & $6.44809E-07$ & 2.6543E-07 & 2.36975E-08 \\ \hline & &Maclaurin on II & Maclaurin on III & Maclaurin on IV \\ \hline\hline & Fourth Order & 5.8797E-07 & 5.83925E-07 & 6.42737E-07 \\ \hline & Fifth Order & 2.30942E-08 & 2.29035E-08 & 2.58348E-08 \\ \hline & Sixth Order & 8.87643E-10 & 8.78999E-10 & 1.01714E-09 \\ \hline\hline Time & R-K on I & R-K on II & R-K on III & R-K on IV \\ \hline\hline $1.25$ & 7.57602E-08 & 1.54213E-07 & 4.93023E-07 & 3.19551E-08 \\ \hline & & Maclaurin on II & Maclaurin on III & Maclaurin on IV \\ \hline\hline & Fourth Order & 1.95583E-08 & 1.51518E-06 & 8.53911E-07 \\ \hline & Fifth Order & 1.03188E-08 & 7.762E-08 & 3.02237E-08 \\ \hline & Sixth Order & 7.6610E-10 & 3.61642E-09 & 1.06356E-09 \\ \hline \end{tabular} \end{center} It is interesting to note that the fourth and fifth order Maclaurin polynomials are similar to the fourth order Runge-Kutta results, but that the sixth order Maclaurin polynomials give noticeable improvement without much labor. It is also interesting to note that the fourth order Runge-Kutta method performs best on (IV) while the same is not true for the Maclaurin polynomials. One of the reasons is that different Maclaurin polynomials are used to update the Maclaurin polynomial for $u$ in \eqref{II}, \eqref{III} and \eqref{IV}. In a future paper we analyze which polynomial projections give the 'best' numerical results. In the next example we show a projection that gives an interesting way to obtain an analytic solution and then we use that system to generate numerical results. \noindent{\bf Example 3.} Consider (I): \begin{equation} \label{3.I} x' = \sin x ; \ x(0) = \alpha. \end{equation} We let $x_2=\sin x$ and $x_3=\cos x$ to obtain (II): \begin{equation} \label{3.II} \begin{gathered} x' = x_2 ; \quad x(0) = \alpha \\ x_2' = x_2x_3; \quad x_2(0) = \sin(\alpha) \\ x_3' = -x_2^2; \quad x_3(0) = \cos(\alpha). \end{gathered} \end{equation} Note that in this system one does not have to use the Maclaurin polynomial for the sine function and that the above algorithm will generate the coefficients of the Maclaurin polynomial for $x$ without knowing the one for sine. We give numerical results for these two systems using fourth order Runge-Kutta and fourth, fifth and sixth order Maclaurin polynomials, but first we show that using the polynomial system it is easy to generate the solution to the original initial value ordinary differential equation. Using the equation for $x_2'$ and $x_3'$ it is seen that $x_2^2+x_3^2=1$ so that $x_3'=x_3^2-1$ whose solution is $$ x_3 = \frac{1-e^{2t+2B}}{1+e^{2t+2B}}. $$ >From this we obtain $$ x_2 = \frac{2e^{t+B}}{1+e^{2t+2B}} $$ and since $ x' = x_2 $, we finally have $$ x = 2 \arctan e^{t+B} -2 \arctan B + \alpha, $$ where $B = \sqrt\frac{1-\cos(\alpha)}{1+\cos(\alpha)}$. We use this exact solution to compare to our numerical results. In the table below is the error of the results of a simulation with a time step of 0.0625 using fourth order Runge-Kutta on the systems \eqref{3.I}--\eqref{3.II} and the fifth order Maclaurin polynomial on \eqref{3.II} using the algorithm above for the initial condition $x(0)=\frac{31\pi}{32}$. The computing times were essentially the same. The errors are in comparison with the exact solution given by Maple using 30 digits of accuracy. \begin{center} \begin{tabular}{|c|c|c|c|c|c|} \hline Time & Order & R-K on I & R-K on II & Maclaurin on II \\ \hline\hline $0.125$ & $4$ & $1.443658E-09$ & 1.35448E-09 & 1.21177E-09 \\ \hline & $5$ & & & 7.6911E-12 \\ \hline & $6$ & & & 6.14E-14 \\ \hline $2$ & $4$ & 3.6696E-09 & 3.56877E-09 & 4.52047E-09 \\ \hline & $5$ & & & 6.74697E-11 \\ \hline & $6$ & & & 1.1585E-12 \\ \hline \end{tabular} \end{center} In these results it is seen that the polynomial system gives the best results and that increasing the degree improves the results significantly. It is well documented that (rational polynomial) Pade' approximants are valuable in approximating analytic functions. Pade' approximants are particularly robust when approximating a wide class of singular functions. However, very little is found in the literature regarding the use of Pade approximants in the numerical solution of differential systems. This is especially true regarding algorithms that operate in a non-symbolic environment. One contributing factor to this is likely the upfront cost of generating the Taylor coefficients of the solution of the system, which are required to create the Pade' approximants. For a given degree, the Taylor coefficients would typically be generated by symbolically differentiating the right hand side of the system multiple times, which characteristically results in computationally complex expressions that need to be evaluated numerous times throughout the solution procedure. Also, if higher degree approximations are desired, the numerical algorithm must be recoded to include the additional expressions first derived symbolically from computing higher order derivatives of the right hand side. As shown in the examples above and documented in \cite{ps1} and \cite{ps2}, the Algebraic-Maclaurin algorithm generates the Taylor coefficients of the solution of the system in a remarkably computationally friendly manner. With this method, generating higher (or lower) degree approximations requires nothing more than a reassignment of a parameter in a running program, which as mentioned above, gives the user ``on the fly" control over the order of the numerical algorithm. The Algebraic-Maclaurin method serves as an excellent vehicle to numerically generate Pade' approximants to the solution of differential systems, and we refer to this method as the Algebraic-Maclaurin-Pade (AMP) method. Consider the following second order nonlinear initial value problem, which can occur in physical examples. A similar equation arises in the {\it N}-body problem, which has been recently studied in \cite{prl} using the results of \cite{ps1} and \cite{ps2}. \noindent{\bf Example 4.} (I): \begin{equation} \label{4.I} x''(t) = x^{3/2} \quad x(0) = 1, \quad x'(0) = 2/ \sqrt{5}. \end{equation} For this equation, the explicit solution is $$ x(t) = \frac{1}{(1-\frac{1}{2\sqrt{5}}t)^4} . $$ Note the singularity in the solution at $t = 2\sqrt{5}$. On letting $x_1=x$, $x_2=x'$ $x_3=x^{1/2}$, and $x_4=\frac{1}{x_3}$, we obtain (II): \begin{equation} \label{4.II} \begin{gathered} x_1' = x_2; \quad x_1(0)= 1 \\ x_2' = x_3^3; \quad x_2(0) = 2/ \sqrt{5} \\ x_3' = \frac{1}{2}x_4x_2;\quad x_3(0)= 1 \\ x_4' = -\frac{1}{2}x_4^3x_2; \quad x_4(0)= 1. \end{gathered} \end{equation} We compare numerical solutions of \eqref{4.I} and \eqref{4.II} using the fourth order Runge-Kutta scheme on \eqref{4.I}, and both the above outlined Algebraic-Maclaurin and AMP methods to generate various order Maclaurin and Rational polynomial approximations for (II). The fourth order Runge-Kutta approximation for (I) and the Maclaurin and Rational polynomial approximations of the first component of (II) are contrasted against the explicit solution. Of course, the Algebraic-Maclaurin approximation matches results from employing a standard Taylor/Maclaurin series numerical method. We present results for relative error and approximate CPU cost for $t$ ``close" to the singularity. If $t=4.4721$ the exact solution gives $x(4.4721)= 2.393\dots \; E+20$. The table below gives the relative error at several orders and number of steps on the interval $[0,4.4721]$. \begin{center} \begin{tabular}{|c||c|c|c|c|c|c|c|} \hline Method & Order & Steps & Relative Error & CPU ($secs$) \\ \hline\hline R-K & 4 & $1.E+05$ & 5.4751E-01 & 3.00E-01 \\ \hline R-K & 4 & $1.E+07$ & 2.0309E-07 & 3.02E+01 \\ \hline R-K & 4 & $1.E+08$ & 2.3338E-08 & 3.03E+02 \\ \hline Algebraic-Maclaurin & 4 & $1.E+05$ & 6.4479E-01 & 1.15E+00 \\ \hline Algebraic-Maclaurin & 32 & $1.E+05$ & 2.4705E-06 & 3.58E+01 \\ \hline Algebraic-Maclaurin & 64 & $1.E+05$ & 3.7209E-09 & 1.36E+02 \\ \hline AMP & 16 & $1.E+03$ & 2.2709E-08 & 2.80E-01 \\ \hline AMP & 36 & $1.E+03$ & 4.1217E-10 & 1.51E+00 \\ \hline \end{tabular} \end{center} Here R-K represents the fourth order Runge-Kutta method, and AMP represents results from the Algebraic-Maclaurin-Pade method. Notice the only parameter that can be changed in the fourth order Runge-Kutta method is the number of steps, but both the Maclaurin, generated by the Algebraic-Maclaurin algorithm, and the AMP method have parameters for both the order and the number of steps. In the results above we have kept the number of steps constant merely to emphasize the impact of the change of order. The results indicate that even though the AMP algorithm is computationally more expensive per step, the increased accuracy, with a considerably longer step size, results in significant improvement in overall $CPU$ costs. It is worth noting that the AMP results above are from an algorithm that uses full Gaussian elimination with scaled partial pivoting to generate Pade' coefficients from known Maclaurin coefficients. There will be a further significant reduction in $CPU$ costs when a numerically stable algorithm which uses a continued fraction technique is developed and replaces the less efficient full Gaussian elimination portion of the AMP algorithm. We could also have posed this IVODE in $P_{9,2}$ by making the projection $x_5=x_3^2,x_6=x_4^2,x_7=x_2x_4,x_8=x_2x_6$ and $x_9=x_3x_5$. \noindent{\bf Example 5.} The proof of Theorem 3 provides an algorithm for ``de-coupling'' any system of differential equations that may be recast as a polynomial system. From Theorem 5 and previous examples, this includes a very wide class of systems of differential equations constructed from elementary functions. E.g. consider the system in $x$ and $y$ \begin{gather*} x' = y+e^{xy} \\ y' = x+2e^{xy}. \end{gather*} To rewrite this in polynomial form, we introduce a new variable $u$, identifying $u$ with $e^{xy}$. The system then becomes \begin{gather*} x' = y+u \\ y' = x+2u \\ u' = u(x^2+2xu+y^2+yu) \end{gather*} We compute $x''$ and $x'''$ in terms of $x, y,$ and $u$ as in the proof of Theorem 3, and write corresponding polynomials $x^{(j)}-P^*_j(x,y,u)$ for $j=0,1,2,3$. Using the polynomial computer algebra system {\sc Singular} (\cite{singular}), we apply elimination techniques to obtain the appropriate Gr\"{o}bner basis polynomial for the following single (polynomial) differential equation in $x$ alone: \begin{align*} &4x^7{x'}+x^6{x'}^2+4x^5{x'}^3+2x^4{x'}^4-4x^3{x'}^5+x^2{x'}^6-4x{x'}^7 +5x^7{x''} \\ &-4x^6{x'}{x''}+9x^5{x'}^2{x''}-4x^4{x'}^3{x''}+3x^3{x'}^4{x''}+4x^2{x'}^5{x''} -x{x'}^6{x''}\\ &+4{x'}^7{x''}-5x^6{x''}^2-10x^4{x'}^2{x''}^2-5x^2{x'}^4{x''}^2-20x^5{x'} +16x^4{x'}^2\\ &+26x^3{x'}^3-12x^2{x'}^4-18x{x'}^5 +4{x'}^6-20x^5{x''}+69x^4{x'}{x''} +4x^3{x'}^2{x''}\\ &-74x^2{x'}^3{x''}+24x{x'}^4{x''}+17{x'}^5{x''}+60x^4{x''}^2-72x^3{x'}{x''}^2 -16x^2{x'}^2{x''}^2\\ &+40x{x'}^3{x''}^2-12{x'}^4{x''}^2-40x^3{x''}^3+18x^2{x'}{x''}^3-8x{x'}^2{x''}^3 +2{x'}^3{x''}^3\\ &+4x^5{x'''}+6x^3{x'}^2{x'''}-4x^2{x'}^3{x'''}+2x{x'}^4{x'''} -4{x'}^5{x'''}+x^4{x''}{x'''}\\ &+4x^3{x'}{x''}{x'''}+4x{x'}^3{x''}{x'''}-{x'}^4{x''}{x'''}+32x^3{x'} -28x^2{x'}^2-20x{x'}^3\\ &+17{x'}^4+20x^3{x''}-48x^2{x'}{x''}+18x{x'}^2{x''}+4{x'}^3{x''}-19x^2{x''}^2 +4x{x'}{x''}^2\\ &+14{x'}^2{x''}^2-2x{x''}^3-8{x'}{x''}^3+{x''}^4-16x^3{x'''}+24x^2{x'}{x'''} +8x{x'}^2{x'''}\\ &-18{x'}^3{x'''}+32x^2{x''}{x'''}-30x{x'}{x''}{x'''} +12{x'}^2{x''}{x'''}+4x{x''}^2{x'''}\\ &-2{x'}{x''}^2{x'''}+4x^2{x'''}^2-4x{x'}{x'''}^2+{x'}^2{x'''}^2-16x{x'} +16{x'}^2-4{x'}{x''}\\ &+16x{x'''}-16{x'}{x'''} +4{x''}{x'''}=0 \end{align*} In an upcoming paper we consider error estimates, convergence rates, and other implementation issues. \subsection*{Acknowledgments} The authors want to thank Professor Darin Stephenson for his helpful conversations \begin{thebibliography}{9} \bibitem{ap} Apostol, Tom M. {\em Calculating Higher Derivatives of Inverses}, The American Mathematical Monthly, 107 : 8 (2000), 738-741. \bibitem{br} Birkhoff, Garrett and Rota, Gian-Carlo. {\em Ordinary Differential Equations, 2nd ed.}, (1969) Blaisdell Publishing Co., Waltham, MA. \bibitem{clo} Cox, David; Little, John; and O'Shea, Donal. {\em Ideals, Varieties, and Algorithms, 2nd ed.} (1996) Springer-Verlag, New York. \bibitem{gw} Gerald, Curtis F. and Wheatley, Patrick O. {\em Applied Numerical Analysis, 6th ed.}, (1999) Addison-Wesley-Longman, Inc., Reading, MA. \bibitem{singular} Greuel, G.-M.; Pfister, G.; and Sch\"onemann, H. {\sc Singular} 2.0. {\em A Computer Algebra System for Polynomial Computations}. Centre for Computer Algebra, University of Kaiserslautern (2001). {\tt http://www.singular.uni-kl.de}. \bibitem{ke} Kendig, Keith. {\em Elementary Algebraic Geometry} (1977) Springer-Verlag, New York. \bibitem{ps1} Parker, G. Edgar, and Sochacki, James S. {\em Implementing the Picard Iteration}, Neural, Parallel, and Scientific Computation, 4 (1996), 97-112. \bibitem{ps2} Parker, G. Edgar, and Sochacki, James S. {\em A Picard-McLaurin Theorem for Initial Value PDE's}, Abstract and Applied Analysis 5 : 1 (2000), 47-63. \bibitem{prl} Pruett, C.D.; Rudmin, J.W.; and Lacy, J.M. {\em An adaptive N-body algorithm of optimal order}, Journal of Computational Physics, 187 (2003), pp. 298-317. \end{thebibliography} \end{document}