\documentclass[11pt]{article}
\usepackage{graphicx,amsmath,amssymb,url}
\usepackage{full page}
\newcommand{\ds}{\displaystyle}

\begin{document}

\title{The Power Series Method}

\author{\textit{compiled by} Stephen Lucas \\ Department of Mathematics and Statistics \\ James Madison University \\ Harrisonburg VA 22807}

\maketitle
%\tableofcontents

\section{Introduction and Motivation} 

A technique introduced in most first courses in differential equations is the Power Series Method (PSM) for linear second order ordinary differential equations with non-constant coefficients. The idea is to substitute a power series for the unknown function and equate coefficients. As long as the differential equation's coefficients are relatively simple, explicit recurrence relations for the power series coefficients can be found and a solution obtained. However, even slightly complicated differential equation coefficients quickly lead to an algebraic morass, and most books quickly move on to a discussion of regular singular points, the method of Frobenius, Bessel functions, and so on. Nonlinear differential equations are avoided, and the method rarely makes a re-appearance. 

A relative of the PSM, known as the class of Taylor methods, is introduced in the field of numerical analysis for approximating the solution of differential equations. A power series is substituted into the differential equation, and this time coefficients of the power series are obtained as functions by successive differentiation of the original differential equation, with extensive application of the chain rule. Typically, the process is truncated at some point, and the first terms of the Taylor series of the solution are obtained. While this approach is straightforward, can deal with nonlinear differential equations, and can be applied on successive intervals leading to a solution on an extended domain, typically the successive differentiations lead to complicated functions and require substantial algebraic manipulation. As a result, they are typically replaced by Runge-Kutta methods that use multiple evaluations of the differential equation functions to match the Taylor series of the solution to a given number of terms. Unfortunately, varying the order of accuracy for Runge-Kutta methods usually requires completely new function evaluations, and finding the solution between data points is far from obvious.

In 1986, Parker and Sochacki [] discovered that if Picard's method is applied to a system of first order differential equations $\vec y' = \vec F(\vec y)$ whose right hand side is polynomial in the dependent variables, then the power series of the solutions emerge, one term per iteration. We shall call this a \textit{polynomial differential system}, and any function that is a solution of a polynomial differential system is \textit{projectively polynomial}. Later, it was shown that the Taylor series coefficients for projectively  polynomial functions can be more easily found by direct substitution using Cauchy products to deal with products of variables -- the PSM. In addition, since the system is autonomous, a truncated Taylor series can be used to approximate the solutions at a later time step, and the process repeated giving solutions over arbitrary intervals.

While it may appear that the applications of the PSM are limited, in fact:
\begin{itemize}
\item Virtually every function analytic at zero is projectively polynomial.
\item Virtually every system of differential equations can be rewritten as a polynomial differential system.
\item Every polynomial differential system can be rewritten in quadratic form, by introducing new variables and either extending the original differential system, or (usually more efficiently) defining the Taylor series coefficients for the new variables in terms of previously calculated coefficients for the original variables.
\item The algorithm to find the Taylor series coefficients for a quadratic polynomial differential system is simple and fast, and as many Taylor coefficients as necessary can be found. The coefficients at the $(n+1)$st order explicitly depend on the previous coefficients, and so extending the order of accuracy is straightforward. Finding a one hundredth order Taylor approximation with fifty digit accuracy is just as easy as a fourth order method.
\item Every quadratic polynomial differential system can be rewritten so that all the quadratic terms are squares of individual variables, with no cross-product terms.
\item Every polynomial differential system can be ``diagonalized'' into differential equations that are polynomial functions of only one variable and its derivatives.
\item There is an a-priori error bound for the solutions of polynomial differential systems, so a combination of step size and order can be chosen to guarantee a given error requirement.
\item The PSM can be easily applied as an ``effectively symplectic'' solver, that conserves some Hamiltonian quantity to machine accuracy. The PSM can also be easily extended to delay differential equations, integro-differential equations, stochastic differential systems, and boundary value problems. It can also be made parallel. It can also be extended to partial differential equations.
\end{itemize} 

Our aim here is to justify these claims, and show how the PSM is a powerful, flexible and efficient approach to solving differential equations. We will also show how the few analytic functions that can't be written as solutions to polynomial differential systems can still be written in a form that is amenable to power series substitution. The level of mathematical sophistication won't rise beyond that of an upper level undergraduate. We will begin by looking at the disadvantages of standard numerical approaches to systems of differential equations, motivating the PSM. We follow this with a detailed description of how to apply the PSM, both theoretical aspects and practical requirements. 

As the author, this is mostly not my own work. I'm just gathering together the work of many, particularly a large number of JMU faculty and students over many years.

\subsection{Linear Differential Equations}

Most first courses on ordinary differential equations include the power series method for solving linear differential equations with non-constant coefficients. Consider the second order differential equation
\[
p_2(x)y''(x)+p_1(x)y'(x)+p_0(x)y(x)=g(x),
\]
with $y(x_0)=c$ and $y'(x_0)=d$, and for now assume that that $p_2(x_0)\ne 0$, so $x_0$ is what is known as an ordinary point. As long as $p_1/p_2$ and $p_0/p_2$ are both analytic at $x_0$, then the differential equation has a unique solution in the neighborhood of $x_0$. To find a solution, assume
\[
y=\sum_{n=0}^\infty a_n(x-x_0)^n,
\]
as well as
\[
y'=\sum_{n=1}^\infty na_n(x-x_0)^{n-1} \quad\mbox{and}\quad y''=\sum_{n=2}^\infty n(n-1) a_n(x-x_0)^{n-2}.
\]
We can substitute these series into the differential equation, and since it is linear, solving for the $a_n$'s lead to linear equations where each $a_n$ can be written in terms of the previous $a_0,a_1,\ldots,a_{n-1}$, and the solution can either be left with $a_0$ and $a_1$ as unknowns, or the initial conditions applied.

For example, consider the very simple $y''+y=0$ with $x_0=0$. Then 
\[
\sum_{n=2}^\infty n(n-1)a_nx^{n-2} + \sum_{n=0}^\infty a_nx^n=0,
\]
or after shifting the index of the second sum,
\[
\sum_{n=2}^\infty \left[ n(n-1)a_n+a_{n-2} \right] x^{n-2}=0.
\]
The right hand side means every power of $x$ on the left has zero coefficient, so
\[
a_n=-\frac{a_{n-2}}{n(n-1)} \quad\mbox{for}\quad n\ge 2. \label{eq:eg1}
\]
As often happens, there is an explicit solution to this recurrence relation. Successive applications give
\begin{align*} a_n & = -\frac{a_{n-2}}{n(n-1)} \ = \ +\frac{a_{n-4}}{n(n-1)(n-2)(n-3)} \\
& = -\frac{a_{n-6}}{n(n-1)(n-2)(n-3)(n-4)(n-5)} \\
& = \cdots = (-1)^{k} \frac{a_{n-2k}}{n(n-1)\cdots(n-2k+2)(n-2k+1)}.
\end{align*}
If $n$ is even, setting $k=n/2$ gives $a_n=(-1)^{n/2}a_0/n!$, and if $n$ is odd, setting $k=(n-1)/2$ gives $a_n=(-1)^{(n-1)/2} a_1/n!$. Combined,
\[
y=a_0\left( 1-\frac{x^2}{2!}+\frac{x^4}{4!}+\frac{x^6}{6!}-\cdots \right) + a_1\left( x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots \right),
\]
which we should recognize as the familiar $y=a_0\cos x+a_1\sin x$.

As a slightly more complicated example, consider $(1-x^2)y''-2xy'+k(k+1)y=0$. Then
\[
(1-x^2) \sum_{n=2}^\infty n(n-1)a_nx^{n-2} - 2x\sum_{n=1}^\infty na_nx^{n-1} + k(k+1)\sum_{n=0}^\infty a_nx^n =0,
\]
or with $m=n-2$,
\[
\sum_{m=0}^\infty (m+2)(m+1)a_{m+2}x^m - \sum_{n=2}^\infty n(n-1)a_nx^n - 2\sum_{n=1}^\infty na_nx^n + k(k+1)\sum_{n=0}^\infty a_nx^n =0,
\]
or
\begin{multline*}
\left( 2\cdot 1 a_2 + k(k+1)a_0 \right)x^0 + \left( 3\cdot 2 a_3 - 2\cdot 1 a_1 + k(k+1)a_1 \right) x^1 \\
+ \sum_{n=2}^\infty \left( (n+2)(n+1)a_{n+2} - (n(n-1)+2n-k(k+1)) a_n \right) x^n =0. \end{multline*}
Equating coefficients with $a_0$ and $a_1$ given, 
\begin{align*}
a_2 & = -\frac{k(k+1)}{2} a_0, \\
a_3 & = \frac{2-k(k+1)}{6} a_1, \quad \mbox{and} \\
a_{n+2} & = \frac{n(n+1)-k(k+1)}{(n+2)(n+1)}a_n, \quad\mbox{for} \quad n=2,3,\ldots.
\end{align*}
If $k$ is a natural number, this last result gives $a_{k+2}=0$, and as long as the initial condition of opposite parity (odd versus even) is chosen zero, the solution is a polynomial. In fact, the classical Legendre polynomials are scaled polynomial solutions of this differential equation.

Some more adventurous students consider problems with more complicated functions as coefficients multiplying the derivatives of $y$, but unless the series combine particularly nicely, the relationships between coefficients can become extremely difficult to work with. In any event, differential equations textbooks quickly move on at this point to Taylor series around \textit{regular singular points}, where $p_2(x_0)=0$ but a unique solution still exists. The \textit{method of Frobenius} finds solutions as power series multiplied by $(x-x_0)^c$ where $c$ is usually not a natural number, but we shall not pursue such problems here. The key concept is that the power series method can be used to solve linear differential equations, particularly with simple coefficients. However, in standard differential equations texts it is never applied to nonlinear equations.

\subsection{Numerical Taylor Methods}

Since solutions of most (systems of) ordinary differential equations can't be found in closed form, an extensive collection of numerical approximation techniques have been developed over the years. Almost every text that mentions the topic starts with \textit{Euler's method}. Given the first order system $\vec y'(t)=\vec f(t,\vec y(t))$ with $\vec y(t_0)=\vec y_0$, approximate the derivatives by $y_i'(t) \approx (y_i(t+\Delta t)-y_i(t))/\Delta t$, then
\[
\vec y(t_0+\Delta t) \approx \vec y_0 + \Delta t f(t_0,\vec y(t_0)).
\]
If we let $t_1=t_0+\Delta t$, then we have approximated the solution at the next time step and can repeat the formula to find approximations to $\vec y(t_1+\Delta t) = \vec y(t_2)$, $\vec y(t_2+\Delta t) = \vec y(t_3)$, and so on, as $\vec y(t_{i+1}) \approx y(t_i) + \Delta t f(t_i,\vec y(t_i))$.

Euler's method is the simplest of a family of methods known as Taylor methods, where the aim is to find the Taylor series of the elements of $\vec y$ at $t=t_i$ to as high an order as possible, then use this Taylor polynomial to approximate $\vec y$ at $t_i+\Delta t=t_{i+1}$ and repeat. Assuming each $y_j(t)$ is infinitely differentiable at $t=t_i$ for $j=1,2,\ldots,n$, then their Taylor series are
\begin{align*}
y_j(t) = &\ y_j(t_i) + \frac{y_j'(t_i)}{1!}(t-t_i) + \frac{y_j''(t_i)}{2!}(t-t_i)^2 + \frac{y_j^{(3)}(t_i)}{3!}(t-t_i)^3 + \cdots + \\
& \frac{y_j^{(n)}}{n!}(t-t_i)^n + \cdots.
\end{align*}
Since $y_j'(t_i)=f_j(t_i,\vec y(t_i))$, and with $t_{i+1}=t_i+\Delta t$, 
\begin{align*}
y_j(t_{i+1}) = &\ y_j(t_i)+f_i(t_i,\vec y(t_i))\frac{\Delta t}{1!} + f_j'(t_i,y(t_i))\frac{\Delta t^2}{2!} +\\
 & f_j''(t_i,y(t_i))\frac{\Delta t^3}{3!} + \cdots + f_j^{(n-1)}(t_i,y(t_i))\frac{\Delta t^n}{n!} + \cdots.
\end{align*}
The Taylor method of order $n$ truncates the Taylor series after the term involving $\Delta t^n$, and the error at each step is of order $\Delta t^{n+1}$. However, when approximating the solution to differential equations on some interval $[a,b]$, the error is considered order $n$, because reducing the size of the time step also increases the number of steps. Thus, if you halve the step size you would expect the error at the end to decrease by a factor of $2^n$.

For example, consider $y'=1+ty$. Remembering that $y$ is a function of $t$ and using the chain rule,
\begin{align*}
f(t,y) & = 1+ty, \\
f'(t,y) & = \frac{d}{dt}f(t,y)=0+y+ty' = y+t(1+ty) = t+(1+t^2)y, \\
f''(t,y) & = \frac{d}{dt}f'(t,y)= 1 + 2ty + (1+t^2)y' \\
 & = 1+2ty + (1+t^2)(1+ty) = 2+t^2+(3t+t^3)y, \mbox{ and} \\
f'''(t,y) & = \frac{d}{dt}f''(t,y) = 2t + (3+3t^2)y + (3t+t^3)y' \\
 & = 2t + (3+3t^2)y + (3t+t^3)(1+ty) = 5t + t^3 + (3+6t^2+t^4)y.
\end{align*}
Thus the fourth order Taylor method for this differential equation is
\begin{align*}
y(t_{i+1}) & \approx y(t_i) + (1+t_iy(t_i))\Delta t + (t_i+(1+t_i^2)y(t_i))\frac{\Delta t^2}{2}  +  \\
 & \quad (2+t_i^2+(3t_i+t_i^3)y(t_i)) \frac{\Delta t^3}{6} + (5t_i + t_i^3 + (3+6t_i^2+t_i^4)y(t_i))\frac{\Delta t^4}{24}.
\end{align*}
The difficulty with Taylor methods is that high order derivatives of $f$ usually become extraordinarily complicated. Even for this simple example, $f'''$ reasonably complicated. As another example, consider the marginally more complicated $y'=1+t^2y^2$. Then
\begin{align*}
f(t,y) & = 1+t^2y^2, \\
f'(t,y) & = 2ty^2+2t^2y+2t^4y^3, \\
f''(t,y) & = 2t^2+8ty+2y^2+8t^4y^2+12t^3y^3+6t^6y^4, \mbox{ and} \\
f'''(t,y) & = 12t+12y+16t^4y+40t^3y^2 + 40t^2y^3+36t^3y^3 + 40t^6y^3+ \\
& \quad 36 t^5y^4 + 36t^5y^5+24t^8y^5.
\end{align*}
Already, the amount of algebra required to find a fourth order Taylor series is substantial and tedious. Another example is $y''=\sin y$, which can be written as the pair of first order equations $y_1'=y_2$ and $y_2'=\sin y_1$. Then
\begin{align*}
f_1& =y_2, & f_2& =\sin y_1, \\
f_1'& =f_2, & f_2'& =y_2\cos y_1, \\
f_1''& =f_1', & f_2''& = -y_2^2\sin y_1+\cos y_1\sin y_1, \\
f_1^{(3)}& =f_2'', & f_2^{(3)}& = -y_2^3\cos y_1-3y_2\sin^2 y_1 + y_2\cos^2 y_1, \\
f_1^{(4)}& =f_2^{(3)}, & f_2^{(4)}& =-y_2^4\sin y_1 -11y_2^2\sin y_1 \cos y_1   \\
& & &  \quad -3\sin^3 y_1  +\cos^2 y_1 \sin y_1.
\end{align*}

Building the Taylor series for even these simple example problems requires substantial effort, and avoiding all of this differentiation is why \textit{Runge-Kutta} methods are one of the most popular approaches to the numerical solution of differential equations. Runge-Kutta methods use multiple function evaluations to form an approximation whose first few terms match those of the actual Taylor series. The disadvantage is that increasing the order of the approximations usually requires a completely new set of function evaluations, and the solution is unavailable between data points. 

\subsection{Power Series Method For These Examples}

The power series method for linear differential equations can be applied to all of these examples, despite the last two not being linear. The first step is to rewrite the differential equations in autonomous form (independent of $t$), which means that a Taylor series can always be expanded around $t=0$: if $\ds \frac{dy}{dt}=f(y(t))$ with $y(a)=y_a$, let $\tau=t-a$. Then $\ds \frac{dy}{dt}=\frac{dy}{d\tau}\frac{d\tau}{dt}$, or $\ds \frac{dy}{d\tau}=f(y(\tau))$ with $y(0)=y_a$. The second step is to make the right hand side polynomial in the independent variables.

Let's reconsider $y'=1+ty$ with initial condition $y(a)=y_a$. By letting $w=t$, we have the first order system $y'=1+wy$ and $w'=1$ with $y(a)=y_a$ and $w(a)=a$. Shifting the time scale by replacing $t$ by $t-a$ gives the equivalent system $y'=1+wy$ and $w'=1$ with $y(0)=y_a$ and $w(0)=a$. Then $w=a+t$, and substituting with $\ds y=\sum_{i=0}^\infty y_i t^i$ into the first differential equation gives 
\[
\sum_{i=0}^\infty (i+1)y_{i+1}t^i = 1+(a+t)\sum_{i=0}^\infty y_i t^i = (1+ay_0) + \sum_{i=1}^\infty (ay_i+y_{i-1}) t^i,
\]
and after equating coefficients (and applying the initial condition) we get that
\[
y_0=y_a, \quad y_1=1+ay_0, \quad\mbox{and} \quad y_{i+1} = \frac{ay_i+y_{i-1}}{i+1} \mbox{ for } i=1,2,\ldots.
\]
Using this recurrence relation to find the Taylor series coefficients for $y$ to arbitrary order is much easier than using the standard Taylor method. And because the system is autonomous, the same recurrence relation can be used to find the Taylor expansion for $y$ at any time, simply replacing $y_a$ by new initial values. More specifically, the coefficients give the Taylor series for the original unshifted problem, truncated as $\ds y(t)\approx \sum_{i=0}^n y_i(t-a)^i$. Letting $t_i=a+\Delta t i$ for some time step $\Delta t$, then $\ds y(t_1) \approx \sum_{i=0}^n y_i \Delta t^i$, the approximation at the next data point. Repeating the recurrence relations with $y_a$ replaced by $y(t_1)$ and implicitly shifted by $t_1$ gives a new set of Taylor coefficients for $\ds y(t)\approx \sum_{i=0}^n y_i(t-t_1)^i$, and $\ds y(t_2)\approx \sum_{i=0}^n y_i\Delta t^i$. Continue for as many steps as necessary. This approach is far easier to use than the original Taylor formula.

Moving on to $y'=1+t^2y^2$ with $y(a)=y_a$, we wouldn't initially consider a series method due to the nonlinearity. However, if we replace $y$ by a power series, then we can form the square making use of the \textit{Cauchy product}, which states that
\[
\left( \sum_{i=0}^\infty a_i t^i \right) \left( \sum_{i=0}^\infty b_i t^i \right) = \sum_{i=0}^\infty \left( \sum_{j=0}^i a_jb_{i-j} \right) t^i. \label{eq:cauchy}
\]
This is easily proven by expanding the left hand side and equating coefficients, and is applicable within the smallest radius of convergence of any of the three power series. As a result, shifting the time scale by replacing $t$ by $t-a$ gives $y'=1+(t+a)^2y^2$ with $y(0)=y_a$, and substituting $\ds y=\sum_{i=0}^\infty y_it^i$ gives
\begin{align*}
\sum_{i=0}^\infty (i+1)y_{i+1}t^i & = 1 + (a^2+2at+t^2) \left( \sum_{i=0}^\infty y_i t^i \right)^2 \\
& = 1 + (a^2+2at+t^2)\sum_{i=0}^\infty \left( \sum_{j=0}^i y_jy_{i-j} \right) t^i.
\end{align*}
The initial condition gives $y_0=y_a$, and equating coefficients and using the Cauchy product we get
\begin{align*}
y_1 & = 1+a^2y_0^2, \\
2y_2 & = 2a^2y_0y_1+2ay_0^2, \quad \mbox{and for } i\ge 2, \\
(i+1)y_{i+1} & = a^2\sum_{j=0}^i y_jy_{i-j} + 2a\sum_{j=0}^{i-1} y_jy_{i-1-j} + \sum_{j=0}^{i-2} y_jy_{i-2-j}.
\end{align*} 
These recurrence relations can be easily applied to find the Taylor series to any order, and the same approach as for the previous example can be used to solve the differential equation over an interval. Also, note that the shifted version can be made autonomous by letting $v=t+a$ and $w=(t+a)^2$, so $y'=1+w^2y^2$, $v'=1$, and $w'=2v$ with $y(0)=y_a$, $v(0)=a$, $w(0)=a^2$ -- a polynomial differential system.

Finally, consider $y''=\sin y$ with $y(0)=y_0$, $y'(0)=y_1$. It initially appears that substituting a power series for $y$ won't lead anywhere useful, since we would need the sine of a power series. However, start by rewriting the differential equation as $y_1'=y_2$, $y_2'=\sin y_1$, with $y_1(0)=y_0$, $y_2(0)=y_1$. In addition, let $y_3=\sin y_1$ and $y_4=\cos y_1$. Then $y_3'=y_1'\cos y_1=y_2y_4$, and $y_4'=-y_1'\sin y_1=-y_2y_3$ with $y_3(0)=\sin y_0$, $y_4=\cos y_0$. So we have the (quadratic) polynomial differential system
\[
\begin{array}{l} y_1'=y_2, \\ y_2' =v_3, \\ y_3' = y_2y_4, \\ y_4'=-y_2y_3, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=y_0, \\ y_2(0) =y_1, \\ y_3(0) = \sin y_0,  \\
y_4(0) = \cos y_0. \end{array}
\]
This polynomial differential system is nonlinear, but only involves products on the right hand side. Substituting the power series $\ds y_1=\sum_{i=0}^\infty a_it^i$, $\ds y_2=\sum_{i=0}^\infty b_it^i$, $\ds y_3=\sum_{i=0}^\infty c_it^i$, and $\ds y_4=\sum_{i=0}^\infty d_it^i$, and using the Cauchy product, we get the explicit recurrence relations
\begin{align*}
a_0&=y_0, & b_0&=y_1, & c_0&=\sin(y_0), & d_0=&\cos(y_0), \\
a_{i+1}&=\frac{b_i}{i+1}, & b_{i+1}&=\frac{c_i}{i+1}, & c_{i+1} &= \frac{1}{i+1} \sum_{j=0}^i b_jd_{i-j}, & d_{i+1}&=\frac{1}{i+1}\sum_{j=0}^i b_jc_{i-j}.
\end{align*}
Finding the Taylor series to arbitrary accuracy is straightforward.

\subsection{Picard's Iteration and the Parker Sochacki Method}

Consider the single first order differential equation $y'(t)=f(t,y(t))$, $y(0)=a$. Integrating on $[0,t]$ and applying the initial condition gives the equivalent integral equation $\ds y(t)=a + \int_0^t f(t,y(t))\,dt$. This can be turned into an iteration as 
\[
y_0(t)=a, \quad\mbox{and}\quad y_{n+1}(t) = a+\int_0^t f(t,y_n(t))\,dt \quad\mbox{for}\quad n=0,1,2,\ldots.
\]
Picard's existence theorem states that the initial value problem has a unique solution near $t=0$ as long as $f$ is Lipschitz continuous in $x$ and continuous in $t$. The proof relies on showing that the integral equation iteration, known as Picard iteration, converges. It is often included in introductory texts on differential equations, but usually with the proviso that it is not of practical worth for solving differential equations, since the integrals quickly become intractable. Extending from a single first order differential equation to a system is straightforward. 

For example, consider $y'=\sin y$ with $y(0)=a$. Picard iterates are
\[
\begin{array}{l}
y_0=a, \\
\ds y_1=a+\int_0^t \sin a\,dt = a+t\sin a, \\
\ds y_2=a+\int_0^t \sin(a+t\sin a)\,dt = a+ \frac{\cos a-\cos(a+t\sin a)}{\sin a}, \\
\ds y_3= a+\int_0^t \sin\left( a+ \frac{\cos a-\cos(a+t\sin a)}{\sin a} \right)\, dt,
\end{array}
\] 
and this last integral cannot be simplified. However,  letting $v=\sin y$ and $w=\cos y$, we have the system $y'=v$, $v'=vw$ and $w'=-v^2$ with $y(0)=a$, $v(0)=\sin a$, $w(0)=\cos a$. Picard iteration for this system with $y_0(t)=a$, $v_0(t)=\sin a$, $w_0(t)=\cos a$ is
\[
\begin{array}{ll}
\ds y_{n+1}(t)=a+\int_0^t v_n(t)\,dt, & \ds v_{n+1}(t)=\sin(a) + \int_0^t v_n(t)w_n(t)\,dt, \\ \ds w_{n+1}(t)=\cos(a)-\int_0^t (v_n(t))^2\,dt, & \mbox{for } n=0,1,2,\ldots. \end{array}
\]
The first few iterations give
\[
\begin{array}{l} \ds y_0=a, \\
\ds y_1=a+t\sin a, \\
\ds y_2=a+t\sin a+\frac{\sin a\cos a}{2} t^2, \\
\ds y_3= a+t\sin a+\frac{\sin a\cos a}{2} t^2 + \frac{\sin a\cos^2 a-\sin^3 a}{6} t^3 - \frac{\sin^3 a\cos a}{12} t^4, \\
\ds y_4=a+t\sin a+\frac{\sin a\cos a}{2} t^2 + \frac{\sin a\cos^2 a-\sin^3 a}{6} t^3 - \frac{(6\cos^2 a-5)\sin a\cos a}{24} t^4 - \\ \ds \quad \quad \frac{(16\cos^3 a-3)\sin^3 a}{120} t^5 - \frac{(2\cos^2 a-1)\sin^3 a\cos a}{36} t^6 - \frac{(4\cos^2 a-3) \sin^3 a\cos^3 a}{252} t^7 + \\ \ds \quad \quad \frac{\sin^5 a\cos^3 a}{504} t^8,
\end{array}
\]
as well as similar polynomials for $v$ and $w$. Since the integrals only involve products, this process can be continued indefinitely, and there does appear to be convergence. In fact, at every step an additional term in the power series for $y$ (and $v$ and $w$) is correct. For example, the terms in $y_4$ up to $t^4$ match the power series for $y$, and in general, terms up to $t^n$ in $y_n$ will match the power series.

The fact that Picard iteration can be applied to polynomial differential systems with as many iterations as required was presented in Parker and Sochacki (1996), where they also proved that each iteration provides an additional term in the power series representation of the solution. In addition, they showed that you don't need to keep the incorrect higher order terms at every step. Returning to our $y'=\sin y$ example, while $y_3$ is a quartic polynomial, we can truncate it to a cubic and with no loss of accuracy continue with $y_4$. Substituting cubics into the iteration leads to sixth degree polynomials, but again only the fourth order terms need to be kept.

The original proof of the correctness of the Parker Sochacki method for the polynomial system $\vec y'=\vec F\left(\vec y\right)$ is difficult to follow. An alternative approach starts by observing that since $\vec F$ is polynomial, substituting the power series solution for $\vec y$ means the coefficients of $t^k$ in $\vec y'$ match those of $\vec F \left( \vec y\right)$. Now rewrite Picard iteration as $\vec y'_{n+1} = \vec F \left( \vec y_n \right)$ with $\vec y(0)=\vec y_0$, where the subscript indicates the iteration. Our proof is inductive, and the initial values are the first terms of the power series for $\vec y$. Now, assume that terms up to $t^n$ in the power series of $\vec y_n$ match those of the true solution $\vec y$. Since $\vec F\left( \vec y_n \right)$ is polynomial in the variables, all the terms up to $t^n$ will match those of $\vec F\left( \vec y\right)$, and so will still match the terms of $\vec y'$ up to $t^n$. In particular, the coefficients of $t^n$ are equal, and on the left due to differentiation involve coefficient of $t^{n+1}$ in $\vec y$. As a result, the coefficients of $\vec y_{n+1}$ up to $t^{n+1}$ will match those of $\vec y$, completing our proof.

While many authors have implemented and extended the Parker Sochacki method, our later experience has been that applying the Power Series Method to polynomial differential systems is a more straightforward way of obtaining the Taylor series coefficients of the solution, particularly when the system is quadratic. It doesn't hurt that the two methods share the same acronym of PSM :-) And finally it is worth mentioning that even for quadratic polynomial differential systems, directly applying the Taylor method is much more cumbersome than the PSM.

\section{The Power Series Method}

While formal power series substitution is usually only seen when solving linear differential equations, we have seen a few examples that show that it is equally applicable when working with polynomial differential systems, and is far easier to work with than classical Taylor methods -- while leading to the same solution. We claim that almost all systems of ordinary differential equations can be rewritten as polynomial differential systems, and formal power series substitution used. Here, we shall see how to convert functions, and hence odes, to polynomial form, how we can always reduce polynomial differential systems to quadratic form, and how the amount of effort required can be often be reduced using intermediate variables that are not represented as the solution of a differential equation. Later, we will see that polynomial differential systems can be rewritten as equivalent odes that are independent, and consider the practicalities of estimating the accuracy of the method.

\subsection{Solving a Quadratic Polynomial Differential System}

Before converting a system of differential equations to polynomial form, it is worth explicitly seeing how straightforward it is to evaluate power series coefficients to arbitrary order. If a system of $n$ equations is quadratic in its variables, then it can be written as
\[
y_i' = a_i + \sum_{j=1}^n b_{ij}y_j + \sum_{j=1}^n \sum_{k=j}^n c_{ijk}y_jy_k \quad\mbox{with}\quad y_i(0)=d_i \quad\mbox{for}\quad i=1,2,\ldots,n,
\]
and each of the $a_i$, $b_{ij}$, $c_{ink}$ and $d_i$ coefficients are constants -- typically most of them zero. Substituting $\ds y_i(t)=\sum_{m=0}^\infty y_{im}t^m$ and equating coefficients leads to the explicit recurrence relations  for the coefficients as
\[
y_{i0}=d_i, \quad y_{i1}=a_i+\sum_{j=1}^n b_{ij}y_{j0} + \sum_{j=1}^n\sum_{k=1}^n c_{ijk} y_{j0}y_{k0},
\]
and
\[
y_{i,m+1}=\frac{1}{m+1}\left[ \sum_{j=1}^n b_{ij}y_{jm} + \sum_{j=1}^n\sum_{k=1}^n c_{ijk} \sum_{p=0}^m y_{jp}y_{k,m-p} \right],
\]
for $i=1,2,\ldots,n$ and $m=1,2,\ldots$. The power series coefficients for each dependent variable are explicit functions of the previously calculated coefficients, and can be calculated to any desired order.

In the case of a system of $n_1$ quadratic differential equations with an additional $n_2$ intermediate variables defined as quadratic combinations of variables in the differential system, we have
\[
y_i' = a_i + \sum_{j=1}^{n_1+n_2} b_{ij}y_j + \sum_{j=1}^{n_1+n_2} \sum_{k=j}^{n_1+n_2} c_{ijk}y_jy_k  \quad\mbox{with}\quad y_i(0)=d_i \quad\mbox{for}\quad i=1,2,\ldots,n_1,
\]
and
\[
y_i = a_i + \sum_{j=1}^{n_1} b_{ij}y_j + \sum_{j=1}^{n_1}\sum_{k=1}^{n_1} c_{ijk}y_jy_k  \quad\mbox{for}\quad i=n_1+1,n_1+2,\ldots,n_1+n_2.
\]
As before, substituting power series and equating coefficients leads to explicit recurrence relations for the coefficients in the order
\[
\begin{array}{lll}
y_{i0}=d_i & \mbox{for} & i=1,2,\ldots,n_1, \\
\ds y_{i0}=a_i+\sum_{j=1}^{n_1} b_{ij}y_{j0} + \sum_{j=1}^{n_1}\sum_{k=1}^{n_1} c_{ijk}y_{j0}y_{k0} & \mbox{for} & i=n_1+1,n_1+2,\ldots,n_1+n_2, \\
\ds y_{i1}=a_i+\sum_{j=1}^{n_1+n_2} b_{ij}y_{j0} + \sum_{j=1}^{n_1+n_2} \sum_{k=1}^{n_1+n_2} c_{ijk} y_{j0}y_{k0} & \mbox{for} & i=1,2,\ldots,n_1, \\
\ds y_{i1}=\sum_{j=1}^{n_1} b_{ij}y_{j1} + \sum_{j=1}^{n_1}\sum_{k=1}^{n_1} c_{ijk}(y_{j0}y_{k1}+y_{j1}y_{k0}) & \mbox{for} & i=n_1+1,n_1+2,\ldots, n_1+n_2,
\end{array}
\]
then
\[
\left\{ \begin{array}{l}
\ds y_{i,m+1}=\frac{1}{m+1}\left[ \sum_{j=1}^{n_1+n_2} b_{ij}y_{jm} + \sum_{j=1}^{n_1+n_2} \sum_{k=1}^{n_1+n_2} c_{ijk} \sum_{p=0}^m y_{jp}y_{k,m-p} \right] \mbox{ for } i=1,2,\ldots,n_1, \mbox{ then} \\
\ds y_{i,m+1} = \sum_{j=1}^{n_1} b_{ij}y_{j,m+1} + \sum_{j=1}^{n_1} \sum_{k=1}^{n_1} c_{ijk} \sum_{p=0}^{m+1} y_{jp}y_{k,m-p} \mbox{ for } i=n_1+1,n_1+2,\ldots,n_1+n_2.
\end{array} \right.
\]
for $m=1,2,\ldots$.

\subsection{Converting to Polynomial Form}

Here we shall see how to rewrite elementary functions as a solution of a polynomial differential system. But first, we note that if the functions $f$ and $g$ are projectively polynomial, so are their sum, product and composition. Given $\vec y'=\vec H(\vec y)$ with $y_1=f$ and $y_2=g$, then:
\begin{itemize}
\item If $v_1=f+g$ then $v_1'=f'+g'=h_1(\vec y)+h_2(\vec y)$. The right hand side is polynomial, so $f+g$ is projectively polynomial.
\item If $v_2=fg$ then $v_2'=f'g+fg'=h_1(\vec y)g + f h_2(\vec y) = h_1(\vec y)y_2 + y_1 h_2(\vec y)$. The right hand side is polynomial, so $fg$ is projectively polynomial.
\item If $v_3(t)=f(g(t))$ then $v_3'=f'(g(t))g'(t) = f'(g(t))h_2(\vec y)$. But the process of forming the polynomial system $H$ decomposes $f(t)$, so following the same pattern with $f(g(t))$ will lead to a polynomial system, and $f\circ g$ is projectively polynomial. In practice it is quite straightforward using the chain rule.
\end{itemize}

\subsubsection{Polynomials and Powers}

\begin{itemize}
\item \textit{Linear}: If $y=at+b$ with constants $a$ and $b$,
\[
y'=a, \mbox{ with } y(0)=b.
\]
\item \textit{Quadratic}: If $y=at^2+bt+c$ with constants $a$, $b$ and $c$, then let $y_1=y$ and $y_2=t$ so $y_1'=2at+b=2ay_2+b$ and $y_2'=1$, so
\[
\begin{array}{l} y_1'=2ay_2+b, \\ y_2'=1, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=c, \\ y_2(0)=0. \end{array}
\]
\item \textit{Cubic}: If $y=at^3+bt^2+ct+d$ with constants $a$, $b$, $c$ and $d$, then let $y_1=y$, $y_2=t$ and $y_3=t^2$, so $y_1'=3at^2+2bt+c = 3ay_3+2by_2+c$, $y_2'=1$, and $y_3'=2t=2y_2$, so
\[
\begin{array}{l} y_1'=2by_2+3ay_3+c, \\ y_2'=1, \\ y_3'=2y_2, , \end{array} \mbox{ with } 
\begin{array}{l} y(0)=d, \\ y_2(0)=0, \\ y_3(0)=0. \end{array}
\]
\item \textit{Polynomial}: If $\ds y=\sum_{i=0}^n a_it^i$ then let $y_1=y$ and $y_{i+1}=t^i$ for $i=1,2,\ldots,n-1$, so $\ds y_1'=\sum_{i=1}^n ia_it^{i-1} = a_1+ \sum_{i=2}^nia_iy_i$, $y_2'=1$, and $y_i' = (i-1)t^{i-2}=(i-1)y_{i-1}$ for $i=3,4,\ldots,n$, so
\[
\begin{array}{lll} \ds y_1'=a_1+\sum_{i=2}^n ia_iy_i, & \mbox{with} & y_1(0)=a_0, \\ 
y_2'=1, & & y_1(0)=0, \\ 
y_i'=(i-1)y_{i-1}, & & y_i(0)=0, \quad \mbox{for } i=3,4,\ldots,n.
\end{array}
\]
Alternatively, if individual powers are required rather than a full polynomial, we could use products. For example, if $y=t^6$, let $y_1=t^6$ and $y_2=t$ with intermediate variables $y_3=t^2$ and $y_4=t^4$. Then
\[
\begin{array}{l} y_1'=6t^5, \\ y_2'=1, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=0, \\ y_2(0)=1, \end{array} \quad\mbox{and}\quad 
\begin{array}{l} y_3=y_2^2, \\ y_4=y_3^2. \end{array}
\]
While more compact, the disadvantage of this approach is that introduces products, and even though most of the coefficients are zero, the PSM will require Cauchy product evaluations. As a result, this approach is not recommended.
\item \textit{Reciprocal}: $y=1/t$ doesn't exist at $t=0$, but if $y=1/(t+a)$ with $a\ne 0$, then $y'=-1/(t+a)^2=-y^2$, so 
\[
\begin{array}{lll} y'=-y^2, & \mbox{with} & y(0)=1/a. \end{array}
\]
\item \textit{General Powers:} If $y=(t+a)^\alpha$ for arbitrary $\alpha$ and $a\ne 0$, then $y'=\alpha(t+a)^{\alpha-1}$. Letting $y_1=y$ and $y_2=1/(t+a)$, then $y_1'=\alpha(t+a)^\alpha (t+a)^{-1}=\alpha y_1y_2$ and $y_2'=-y_2^2$, so
\[
\begin{array}{lll} y_1'=\alpha y_1y_2, & \mbox{with} & y_1(0)=a^\alpha, \\
y_2'=-y_2^2, & & y_2(0)=1/a. \end{array}
\]
\item \textit{Reciprocal of Function:} If $f(t)$ is projectively polynomial, it is a solution of a polynomial system of differential equations, which can be written to include $f'=g$ for some function $g$ which is also in the system of differential equations. If we wish to represent $y=1/f$, then by the chain rule, $y'=(-1/f^2)f'=-y^2g$ with $y(0)=1/z(0)$. We then have some choices on how to reduce this to quadratic form.
\begin{itemize}
\item If we let $y_1=y$ and $y_2=y_1g$, then $y_1'=-y_1y_2$ and $y_2'=y_1'g+y_1g' = -y_1y_2g+y_1g' = -y_2^2+y_1g'$. So
\[
\begin{array}{lll} y_1'=-y_1y_2 & \mbox{with} & y_1(0)=1/z(0), \\
y_2'=-y_2^2 + y_1g' & & y_2(0)=g(0)/z(0). \end{array}
\]
This system is quadratic with three products, and so three Cauchy products are required in finding the relevant power series.
\item Allowing the use of intermediate variables, let $y_1=y$ and $y_2=y_1^2$ as before. Then
\[
\begin{array}{lll} y_1'=y_2g & \mbox{with} &  y_1(0)=1/f(0), \end{array} \quad \mbox{and} \quad y_2=y_1^2.
\]
Only two Cauchy products are required to form the power series this way, and is the recommended approach.
\item If we don't insist on a differential system, and already have the series $\ds f(t)=\sum_{i=0}^\infty f_it^i$, then $y=1/f$ can be rewritten as $yf=1$, which involves a single Cauchy product. Then 
\[
\left( \sum_{i=0}^\infty y_i t^i \right) \left( \sum_{i=0}^\infty f_it^i \right) = \sum_{i=0}^\infty \left( \sum_{j=0}^i y_if_{j-i} \right) t^i =1,
\]
so $y_0f_0=1$ and $\ds \sum_{j=0}^i y_if_{j-i}=0$, or
\[
y_0=\frac{1}{f_0} \quad\mbox{and}\quad y_i = -\frac{1}{f_0}\sum_{j=0}^{i-1} y_if_{j-i} \quad\mbox{for } i=1,2,\ldots .
\]
This is the most efficient way to find the Taylor series of the reciprocal of a Taylor series, but has the disadvantage of being a special case.
\end{itemize}

\item \textit{Function Powers:} If $f(t)$ is projectively polynomial, satisfies $f'=g$, and $f(0)\ne 0$, then let $y_1=f^\alpha$ and $y_2=1/f$. Then $y_1'=\alpha f^{\alpha-1}f' = \alpha y_1y_2g$ and $y_2'=(-1/f^2)f'=-y_2^2g'$, with $y_1(0)=f(0)^\alpha$ and $y_2(0)=1/f(0)$. This pair of polynomial differential equations requires four Cauchy products, which can be reduced to three using the intermediate variable $y_3=y_2g$. So
\[
\begin{array}{lll} y_1' = \alpha y_1y_3 & \mbox{with} & y_1(0)=(f(0))^\alpha, \\
y_2'=-y_2y_3 & & y_2(0)=1/f(0), \end{array}
\quad\mbox{and}\quad y_3=y_2g.
\]
We can make the system quadratic using $y_3'=y_2'g+y_2g' = -y_2y_3g+y_2g' = -y_3^2+y_2g'$. So
\[
\begin{array}{lll} y_1'=\alpha y_1y_3, & \mbox{with} & y(0)=(f(0))^\alpha, \\
y_2'=-y_2y_3, & & y_2(0)=1/f(0), \\
y_3'=-y_3^2+y_2g', & & y_3(0)=g(0)/f(0). \end{array}
\]
This version requires four Cauchy products.

As for reciprocal, we can do better if we don't use polynomial projection. Differentiating $y=f^\alpha$ gives $y'=\alpha f^{\alpha-1}f'$, or $y'f=\alpha f^\alpha f'=\alpha yf'$, or
\[
\left( \sum_{i=0}^\infty (i+1)y_{i+1}t^i \right) \left( \sum_{i=0}^\infty f_it^i \right) = \left( \sum_{i=0}^\infty y_it^i \right) \left( \sum_{i=0}^\infty (i+1)f_{i+1}t^i \right).
\]
Expanding both products and equating coefficients of $t^i$ gives
\[
\sum_{j=0}^i (i-j+1)y_{i-j+1}f_j = \sum_{j=0}^i (j+1)y_{i-j}f_{j+1},
\]
or explicitly,
\[
y_0=f_0^\alpha \quad\mbox{and}\quad y_{i+1}=\frac{\alpha y_0f_{i+1}}{f_0} + \frac{(\alpha-1)}{(i+1)f_0}\sum_{j=1}^i jf_jy_{i-j+1} \quad\mbox{for } i=0,1,2,\ldots.
\]
\item \textit{Special Case}: If $y'=y^\alpha$ with $y(0)=\beta\ne 0$, then let $y_1=y$, $y_2=y^\alpha$ and $y_3=1/y$ with $y_2(0)=\beta^\alpha$ and $y_3(0)=1/\beta$, so $y_1'=y^\alpha=y_2$,  $y_2'=\alpha y^{\alpha-1} y' = \alpha y^{2\alpha-1} = \alpha y_2^2y_3$ and $y_3'=(-1/y^2)y'=-y_1^{\alpha-2}=-y_2y_3^2$. This version requires four Cauchy products, but the product $y_2y_3$ is repeated. If we let$y_4=y_2y_3=y^{\alpha-1}$ with $y_4(0)=\beta^{\alpha-1}$, then $y_1'=y_2$, $y_2'=\alpha y_2y_4$, $y_3'=-y_3y_4$ and $y_4'=(\alpha-1)y^{\alpha-2}y' = (\alpha-1)y^{2\alpha-2}=(\alpha-1)y_4^2$. This version requires three Cauchy products, but we can do even better by realizing that we don't actually need $y_3$ in the other equations, and $y_2=y_1y_4$, so with only two Cauchy products,
\[
\begin{array}{lll}
y_1=y & y_1'=y_1y_4, & y(0)=\beta, \\
y_4=y^{\alpha-1}, & y_4'=(\alpha-1)y_4^2, & y_4(0)=\beta^{\alpha-1}.
\end{array}
\]
If, however, $\alpha=(n-1)/n$ for natural number $n>1$, then 
\end{itemize}

\subsubsection{Exponential and Logarithmic}

\begin{itemize}
\item If $y=e^t$, then $y'=e^t=y$ with $y(0)=1$.
\item If $y=e^{at+b}$, then $y'=ae^{at+b}=ay$ with $y(0)=e^b$.
\item If $y=e^{f(t)}$ and $f'=g$, then $y'=e^{f}f'=yg$ with $y(0)=e^{f(0)}$.
\item If $y'=e^y$ then let $y_1=y$ and $y_2=e^y$, so $y_1'=y_2$ and $y_2'=e^yy'=y_2^2$ with $y_1(0)=y(0)$ and $y_2(0)=e^{y(0)}$.
\item If $y=\ln(a+t)$, $a>0$, let $y_1=y$ and $y_2=1/(a+t)$. Then 
\[
\begin{array}{lll} 
y_1'=y_2, & \mbox{with} & y_1(0)=\ln(a), \\
y_2'=-y_2^2, & & y_2(0)=1/a. \end{array}
\]
\item If $y=\ln(f(t))$, $f(0)>0$ and $f'=g$, then let $y_1=y$ and $y_2=1/f(t)$, with the intermediate variable $y_3=y_2g$. Then
\[
\begin{array}{l} y_1'=y_3, \\ y_2'=-y_2y_3, \end{array} \mbox{ with } 
\begin{array}{l} y_1(0)=\ln(f(0)), \\ y_2(0)=1/f(0), \end{array}
\quad\mbox{ and }\quad
y_3=y_2g.
\]
\item If $y=f(t)^{g(t)}$ with $f'=h$ and $g'=i$, then $y'=f^g(g'\ln f+gf'/f)$. Combining various of the results so far, let $y_1=y$, $y_2=\ln(f)$ and $y_3=1/f$. Then $y_1'=y_1y_2i+y_1y_2gh$, $y_2'=y_3h$ and $y_3'=-y_3^2h$ is polynomial. Including the intermediate variables $y_4=y_2i$, $y_5=y_3h$ and $y_6=y_1g$ gives the quadratic system
\[
\begin{array}{l} y_1'=y_1y_4+y_5y_6, \\ y_2'=y_5, \\ y_3=-y_5h, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=f(0)^g(0), \\ y_2(0)=\ln(f(0)), \\ y_3(0) = 1/f(0), \end{array} \quad\mbox{and}\quad 
\begin{array}{l} y_4=y_2i, \\ y_5=y_3h, \\ y_6=y_1g. \end{array}
\]
Finding the coefficients requires 6 Cauchy products.
\end{itemize}

\subsubsection{Trigonometric Functions}

\begin{itemize}
\item To find power series for either $\cos(at+b)$ or $\sin(at+b)$, let $y_1=\cos(at+b)$ and $y_2=\sin(at+b)$ with $y_1(0)=\cos b$ and $y_2(0)=\sin b$. then
\[
\begin{array}{l} y_1'=-ay_2, \\ y_2'=ay_1, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=\cos b, \\ y_2(0)=\sin b. \end{array}
\]
Alternatively, if sine and cosine come in pairs, let $y_1=\alpha\cos(at+b)+\beta\sin(at+b)$ and $y_2=-\alpha\sin(at+b)+\beta\cos(at+b)$, so
\[
\begin{array}{l} y_1'=ay_2', \\ y_2'=-ay_1, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=\alpha\cos b+\beta\sin b, \\ y_2(0)=-\alpha\sin a + \beta\cos b. \end{array}
\] 
\item To find power series for $\cos(f(t))$ or $\sin(f(t))$ given $f'=g$, let $y_1=\cos(f(t))$ and $y_2=\sin(f(t))$. Then
\[
\begin{array}{l} y_1'=-y_2g, \\ y_2'=y_1g, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=\cos(f(0)), \\ y_2(0)=\sin(f(0)). \end{array}
\]
\item Given $y=\tan t$, we have two good alternatives. Since $y'=\sec^2 t$, and $(\sec t)' = \sec t \tan t$,
\[
\begin{array}{l} y_1=\tan t, \\ y_2=\sec t, \end{array} \quad
\begin{array}{l} y_1'=y_2^2, \\ y_2'=y_1y_2, \end{array} \mbox{ with }
\begin{array}{l} y_1(0)=1, \\ y_2(0)=1, \end{array}
\]
with two Cauchy products. This version is useful if the power series for $\sec x$ is also needed in the system. Or, since $\sec^2 t = 1+\tan^2 t$, with only one Cauchy product,
\[
y'=1+y^2, \mbox{ with } y(0)=1.
\]
If $y_1=\tan(f(t))$ with $f'=g$, the use of one intermediate variable gives
\[
y_1'=y_2g, \mbox{ with } y_1(0)=\tan(f(0)), \quad\mbox{and}\quad y_2=1+y_1^2
\]
\item If with $a\ne 0$ we need $\csc(a+t)$, or $\csc(a+t)$ and $\cot(a+t)$, since $(\cot (a+t))'=-\csc^2(a+t)$ and $(\csc(a+t))'=-\csc(a+t) \cot(a+t)$,
\[
\begin{array}{l} y_1=\cot t, \\ y_2=\csc t, \end{array}, \quad
\begin{array}{l} y_1'=-y_2^2, \\ y_2'=-y_1y_2, \end{array} \mbox{ with }
\begin{array}{l} y_1(0)=\cot a, \\ y_2(0)=\csc a. \end{array}
\]
Or, only for $\cot(a+t)$ using $\csc^2(a+t) = 1+\cot^2(a+t)$,
\[
y'=-1-y^2, \mbox{ with } y(0)=\cot a,
\]
\item If $y_1=\sin^{-1}t$, then $y_1'=1/\sqrt{1-t^2}$, and using the earlier polynomial and power results suggests $y_2=1-t^2$, $y_3=t$, $y_4=1/\sqrt{1-t^2}$ and $y_5=1/(1-t^2)$, as well as the intermediate variable $y_6=y_3y_5$,
\[
\begin{array}{l} y_1'=y_4, \\ y_2'=-2y_3, \\ y_3'=1, \\ y_4'=y_4y_6, \\ y_5'=2y_5y_6, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=0, \\ y_2(0)=1, \\ y_3(0)=0, \\ y_4(0)=1, \\ y_5(0)=1, \end{array}
\quad\mbox{and}\quad  y_6=y_3y_5.
\]
This quadratic differential system requires three Cauchy products, and is relatively complicated. Alternatively, given $y_1=\sin^{-1}t$ or $\sin y_1=t$, differentiating with respect to $t$ gives $y_1'\cos y_1=1$ or $y_1'=\sec y$. Letting $y_2=\sec y$ and $y_3=\tan y$, as well as the intermediate variable $y_4=y_2^2$, then
\[
\begin{array}{l} y_1'=y_2, \\ y_2'=y_3y_4, \\ y_3'=y_2y_4, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=0, \\ y_2(0)=1, \\ y_3(0)=0, \end{array}
\quad\mbox{and}\quad y_4=y_2^2.
\]
Still three Cauchy products, but a simpler representation. Also, $\cos^{-1} x = \pi/2 - \sin^{-1}x$, and inverse sine of a function follows as in previous examples.
\item If $y_1=\tan^{-1}t$, let $y_2=\cos y_1$ and $y_3=\sin y_1$, with intermediate variable $y_4=y_2^2$. Then
\[
\begin{array}{l} y_1'=y_4, \\ y_2'=-y_3y_4, \\ y_3'=y_2y_4, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=0, \\ y_2(0)=1, \\ y_3(0)=0, \end{array}
\quad\mbox{and}\quad y_4=y_2^2.
\]
\item If $y_1=\sec^{-1}(t+a)$, let $y_2=\cos y_1$, $y_3=\sin y_1$, $y_4=\cot y_1$ and $y_5=\csc y_1$, with intermediate variables $y_6=y_2y_4$ and $y_7=y_5y_6$. Then
\[
\begin{array}{l} y_1'=y_6, \\ y_2'=-y_3y_6, \\ y_3'=y_2y_6, \\ y_4'=-y_5y_7, \\ y_5'=-y_4y_7, \end{array} \mbox{ with } \begin{array}{l} y_1(0)=a, \\ y_2(0)=\cos a, \\ y_3(0)=\sin a, \\ y_4(0)=\cot a, \\ y_5(0)=\csc a, \end{array}
\quad\mbox{and}\quad \begin{array}{l} y_6=y_2y_4, \\ y_7=y_5y_6. \end{array}
\]
Alas, six Cauchy products for this surprisingly difficult function.
\end{itemize}

\subsubsection{More Complicated Examples}

More complicated compositions or products of functions, either by themselves or as part of a differential system, can be rewritten as quadratic differential systems by successively applying the above rules. One approach is from the outside-in, which  replaces entire functions by new variables, and rewrites their derivatives in terms of previously calculated functions when possible.

For example, consider $y_1=e^t\sin t$ which satisfies $y_1'=e^t\sin t+e^t\cos t=y_1+e^t\cos t$ with $y_1(0)=0$. Letting $y_2=e^t\cos t$ means $y_1'=y_1+y_2$, and $y_2'=e^t\cos t - e^t\sin t = y_2-y_1$ with $y_2(0)=1$. So
\[
\begin{array}{l} y_1=e^t\sin t, \\ y_2=e^t\cos t, \end{array} \mbox{becomes}\quad 
\begin{array}{l} y_1'=y_1+y_2, \\ y_2'=y_2-y_1, \end{array} \mbox{with} \begin{array}{l} y_1(0)=0, \\ y_2(0)=1. \end{array}
\]
The case $y_1'=e^{y_1}\sin y_1$ with $y_1(0)=a$ is only slightly more complicated. Letting $y_2=e^{y_1}\sin y_1$ and $y_3=e^{y_1}\cos y_1$, then $y_1'=y_2$, $y_2'=y_2^2+y_2y_3$ and $y_3'=-y_2^2+y_2y_3$ with $y_1(0)=a$, $y_2(0)=e^a\sin a$ and $y_3(0)=e^a\cos a$. While this system is quadratic, products are duplicated, so to minimize effort add the intermediate variables $y_4=y_2^2$ and $y_5=y_2y_3$. Then
\[
\begin{array}{l} y_1'=e^{y_1}\sin y_1, \\ y_2=e^{y_1}\sin y_1, \\ y_3=e^{y_1}\cos y_1, \end{array} \mbox{becomes} \quad 
\begin{array}{l} y_1'=y_2, \\ y_2'=y_4+y_5, \\ y_3'=-y_4+y_5, \end{array} \mbox{with} \begin{array}{l} y_1(0)=a, \\ y_2(0)=e^a\sin a, \\ y_3(0)=e^a\cos a, \end{array} \quad\mbox{and}\quad
\begin{array}{l} y_4=y_2^2, \\ y_5=y_2y_3. \end{array}
\]

As another example, consider the function $y_1=\sin(\exp(t^2-7))$, which satisfies $y_1'=2t\exp(t^2-7)\cos(\exp(t^2-7))$ with $y_1(0)=\sin(\exp(-7))$. Letting $y_2=2t$, $y_3=\exp(t^2-7)$ and $y_4=\cos(\exp(t^2-7))$, then $y_1'=y_2y_3y_4$, $y_2'=2$, $y_3'=2t\exp(t^2-7) = y_2y_3$, and $y_4'=-2t\exp(t^2-7)\sin(\exp(t^2-7)) = -y_1y_2y_3$ with $y_2(0)=0$, $y_3(0)=\exp(-7)$ and $y_4(0)=\cos(\exp(-7))$. To reduce this system to quadratic, introduce the intermediate variable $y_5=y_2y_3$, so
\[
\begin{array}{l} y_1=\sin(\exp(t^2-7)), \\ y_2=t, \\ y_3=\exp(t^2-7), \\ y_4=\cos(\exp(t^2-7), \end{array} \mbox{becomes}\quad
\begin{array}{l} y_1'=2y_4y_5, \\ y_2'=2, \\ y_3'=y_5, \\ y_4'=-y_1y_5, \end{array} \mbox{with} \begin{array}{l} y_1(0)=\sin(\exp(-7)), \\ y_2(0)=0, \\ y_3(0)=\exp(-7), \\ y_4(0)=\cos(\exp(-7)), \end{array} \quad\mbox{and}\quad y_5=y_2y_3. 
\]

\subsubsection{Cross Terms}

Consider $y_1'=y_1^ay_2^b$ and $y_2'=y_1^cy_2^d$ with arbitrary $a,b,c,d$, and $y_1(0)=\alpha \ne0$, $y_2(0)=\beta\ne 0$. Letting $y_3=y_1^{a-1}y_2^b$, $y_4=y_1^cy_2^{d-1}$ and using the intermediate variable $y_5=y_3y_4$,
\[
\begin{array}{l} y_1'=y_1y_3, \\ y_2'=y_2y_4, \\ y_3'=(a-1)y_3^2+by_5, \\ y_4'=cy_5+(d-1)y_4^2, \end{array} \mbox{with} \begin{array}{l} y_1(0)=\alpha, \\ y_2(0)=\beta, \\ y_3(0)=\alpha^{a-1}\beta^b, \\ y_4(0)=\alpha^c\beta^{d-1}, \end{array}
\quad\mbox{and}\quad y_5=y_3y_4.
\]

\subsection{Special Cases}

Show the problem with $(f(x)-f(x_0))/(x-x_0)$, proofs not polynomial, show how 

\subsection{All Differential Systems are Quadratic}

For example, consider the system 
\begin{equation}
y_1'=y_1+y_2y_3^4, \quad y_2'=y_2y_3, \quad y_3'=y_1^3+y_2^2y_3.
\end{equation}
Let $v_{ijk}=y_1^iy_2^jy_3^k$ for $0\le i\le 3$, $0\le j\le 2$, $0\le k\le 4$ with at least one of $i,j,k>0$. Then there are an additional $4\cdot 3\cdot5-1=59$ valid additional variables, and we want to find a subset of them such that a system of first order odes includes the above three. Then following the proof algorithm \newpage
\[
v_{100}'=v_{100}+v_{014}, \quad v_{010}'=v_{011}, \quad v_{001}'=v_{300}+v_{021},
\]
with
\[
\begin{array}{lll}
v_{014}' & = & v_{004}y_2'+4v_{013}y_3' = v_{004}v_{011} + 4v_{013}(v_{300} + v_{021}) \\
v_{011}' & = & v_{001}y_2'+v_{010}y_3' = v_{001}v_{011} + v_{010}(v_{300} + v_{021}) \\
v_{300}' & = & 3v_{200}y_1' = 3v_{200}(v_{100}+v_{014}) \\
v_{021}' & = & 2v_{011}y_2'+v_{020}y_3' = 2v_{011}v_{011} + v_{020}(v_{300}+v_{021}) \\
v_{004}' & = & 4v_{003}y_3' = 4v_{003}(v_{300}+v_{021}) \\
v_{013}' & = & v_{003}y_2' + 3v_{012}y_3' = v_{003}v_{011} + 3v_{012}(v_{300} +v_{021}) \\
v_{200}' & = & 2v_{100}y_1' = 2v_{100}(v_{100}+v_{014}) \\
v_{020}' & = & 2v_{010}y_2' = 2v_{010}v_{011} \\
v_{003}' & = & 3v_{002}y_3' = 3v_{002}(v_{300}+v_{021}) \\
v_{012}' & = & v_{002}y_2' + 2v_{011}y_3' = v_{002}v_{011} + 2v_{011}(v_{300}+v_{021}) \\
v_{002}' & = & 2v_{001}y_3' = 2v_{001}(v_{300}+v_{021}).
\end{array}
\]
To make this quadratic means a total of $14$ equations, and requires a total of $16$ Cauchy products.

We can cut down on the number of products by looking for duplicates. For example, notice that $v_{013}'$ includes the product $v_{003}v_{011} = (y_3^3)(y_2y_3)=y_2y_3^4=v_{014}$. As a result, we can eliminate the equations for $v_{003}$. Similarly $v_{001}v_{011}=v_{012}$ and $v_{002}v_{011}=v_{013}$,  giving the system
\[
\begin{array}{lll}
v_{100}' & = & v_{100}+v_{014} \\
v_{010}' & = & v_{011} \\
v_{001}' & = & v_{300}+v_{021} \\
v_{014}' & = & v_{004}v_{011} + 4v_{013}(v_{300} + v_{021}) \\
v_{011}' & = &  v_{012} + v_{010}(v_{300} + v_{021}) \\
v_{300}' & = & 3v_{200}(v_{100}+v_{014}) \\
v_{021}' & = & 2v_{011}v_{011} + v_{020}(v_{300}+v_{021}) \\
v_{004}' & = & 4v_{003}(v_{300}+v_{021}) \\
v_{013}' & = & v_{014} + 3v_{012}(v_{300} +v_{021}) \\
v_{200}' & = & 2v_{100}(v_{100}+v_{014}) \\
v_{020}' & = & 2v_{010}v_{011} \\
v_{012}' & = & v_{013} + 2v_{011}(v_{300}+v_{021}) \\
v_{002}' & = & 2v_{001}(v_{300}+v_{021}). \\
\end{array}
\]
This version has $13$ equations and $11$ Cauchy products. It may be possible to reduce the number of Cauchy products further, but there are no obvious approaches to achieve this.

Better: use auxiliary variables. Then we have
\[
y_1'=y_1+v_4y_3, \quad y_2'=v_1, \quad y_3'=v_3y_1+y_2v_1,
\]
with
\[
v_1=y_2y_3, \quad v_2=y_3^2, \quad v_3=y_3^2, \quad v_4=v_1v_2,
\]
and there are only $7$ Cauchy products, the minimum possible.


\subsection{Practicalities}

\subsection{Separating Variables}

\section{A Priori Error Analysis}

\section{Historical Relationships}

\subsection{Fehlberg}

\subsection{AD}

\subsection{Adomium}

\subsection{Others}

\section{Hamiltonian Systems and Symplectic Solvers}

\section{Delay Differential Equations}

\section{Intero-Differential Equations}

\section{High Accuracy Systems and Stiffness}

\section{Specific Systems}

\subsection{N Body Problem}

\subsection{Double Pendulum}

\subsection{Heteroclinic Systems}

\subsection{Stochastic Differential Equations}

\subsection{Predator-Prey}

\subsection{Leah Cosine and Sine}

\section{Other Applications}

\subsection{Newton's Method}

\subsection{Boundary Value Problems}





\begin{thebibliography}{99}
\bibitem{Parker} Parker \& Sochacki
\end{thebibliography}


\end{document}

