-------------------------------------------------
lines 7-108 of file: example/chkpoint_two/ode.cpp
-------------------------------------------------

{xrst_begin chkpoint_two_ode.cpp}
{xrst_spell
   kutta
   runge
}

Checkpointing an ODE Solver: Example and Test
#############################################

Purpose
*******
In this example we :ref:`checkpoint<chkpoint_two-name>` one step of an ODE solver.

Problem
*******
We consider the initial value problem with parameter :math:`x` defined by,
:math:`z(0, x) = z_0 (x)`,

.. math::

   \partial_t z(t, x ) = h [ x , z(t, x) ]

Note that if :math:`t` needs to be in the equation, one can define
the first component of :math:`z(t, x)` to be equal to :math:`t`.

ODE Solver
**********
For this example, we consider the Fourth order Runge-Kutta ODE solver.
Given an approximation solution at time :math:`t_k` denoted by
:math:`\tilde{z}_k (x)`, and :math:`\Delta t = t_{k+1} - t_k`,
it defines the approximation solution :math:`\tilde{z}_{k+1} (x)`
at time :math:`t_{k+1}` by

.. math::
   :nowrap:

   \begin{eqnarray}
   h_1 & =  & h [ x , \tilde{z}_k (x) ]
   \\
   h_2 & =  & h [ x , \tilde{z}_k (x) + \Delta t \; h_1 / 2 ]
   \\
   h_3 & =  & h [ x , \tilde{z}_k (x) + \Delta t \; h_2 / 2 ]
   \\
   h_4 & =  & h [ x , \tilde{z}_k (x) + \Delta t \; h_3 ]
   \\
   \tilde{z}_{k+1} (x) & = &
      \tilde{z}_k (x) + \Delta t \; ( h_1 +  2 h_2 + 2 h_3 + h_4 ) / 6
   \end{eqnarray}

If :math:`\tilde{z}_k (x) = z_k (x)`,
:math:`\tilde{z}_{k+1} (x) = z_{k+1} (x) + O( \Delta t^5 )`.
Other ODE solvers can use a similar method to the one used below.

ODE
***
For this example the ODE is defined by
:math:`z(0, x) = 0` and

.. math::

   h[ x, z(t, x) ] =
   \left( \begin{array}{c}
         x_0                     \\
         x_1 z_0 (t, x)          \\
         \vdots                  \\
         x_{n-1} z_{n-2} (t, x)
   \end{array} \right)
   =
   \left( \begin{array}{c}
         \partial_t z_0 (t , x)      \\
         \partial_t z_1 (t , x)      \\
         \vdots                      \\
         \partial_t z_{n-1} (t , x)
   \end{array} \right)

Solution
********
The solution of the ODE for this example,
which is used to check the results,
can be calculated by
starting with the first row and then using the solution
for the first row to solve the second and so on.
Doing this we obtain

.. math::

   z(t, x) =
   \left( \begin{array}{c}
      x_0 t                  \\
      x_1 x_0 t^2 / 2        \\
      \vdots                 \\
      x_{n-1} x_{n-2} \ldots x_0 t^n / n !
   \end{array} \right)

Source
******
{xrst_literal
   // BEGIN C++
   // END C++
}

{xrst_end chkpoint_two_ode.cpp}