notes.tex

\documentclass[10pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{epstopdf}
\usepackage{abstract}
\usepackage[toc,page]{appendix}
\usepackage[sort]{cite}

\usepackage{gensymb}
\usepackage{hyperref}


\newcommand{\half}[0]{\frac{1}{2}}
\newcommand{\bvec}[1]{\mathbf{#1}}
\newcommand{\bigO}[2]{\mathcal{O}\left(#1^{#2} \right)}
\newcommand{\dotprod}[2]{ \left<#1 , #2 \right> }

\newcommand{\dd}[0]{ \mathrm{d} }

\title{ Some notes on implicit Runge-Kutta methods }
\author{ Stefan Paquay }
\date{  }

\begin{document}
\maketitle

\section{The idea}
Runge-Kutta methods are multi-stage, single-step methods aimed towards solving (systems of) ODEs.
They work by constructing at each time step approximations to the derivative of the function in between the current time level $t$ and the next $t + \Delta t,$ after which a linear combination of these stages is used to advance the numerical approximation to $y$ to the next level.
That is, if we have $\dd y/ \dd t = f(t,y),$ and $y_n$ is the numerical approximation to $y$ at $t_n,$ then we have
\begin{equation*}
  k_i = f\left( t + c_i \Delta t, y_n + \Delta t \sum_{j=0}^{N-1} \alpha_{i,j} k_j \right), \qquad y_{n+1} = y_n + \Delta t \sum_{i=0}^{N-1} b_i k_i.
\end{equation*}
The methods can be summarized easily in so-called Butcher tableaus, which conveniently list $b_i,~c_i$ and $\alpha_{i,j}.$ See \ref{tab:butcher} for the Butcher tableau of some methods.

\begin{table}[b]
  \centering
  \caption{Butcher tableaus for some Runge-Kutta methods: The implicit midpoint (left), the classical Runge-Kutta method (center), and the fourth order Gauss-Legendre method (right). \label{tab:butcher}}
  \def\arraystretch{1.5}%
  \begin{tabular}{c|c}
    $\half$ & $\half$ \\ \hline
    & 1
  \end{tabular} \quad
  \begin{tabular}{c|cccc}
    0 & & & & \\
    $\half$ & $\half$ & & & \\
    $\half$ & & $\half$ & & \\
    1 & & & 1 & \\ \hline
    {} & $\frac{1}{6}$ & $\frac{2}{6}$ & $\frac{2}{6}$ & $\frac{1}{6}$
  \end{tabular}
  \begin{tabular}{c|cc}
    $\half - \frac{\sqrt{3}}{6}$ & $\frac{1}{4}$ & $\frac{1}{4} - \frac{\sqrt{3}}{6}$ \\
    $\half + \frac{\sqrt{3}}{6}$ & $\frac{1}{4} + \frac{\sqrt{3}}{6}$ & $\frac{1}{4}$ \\ \hline
    {} & $\half$ & $\half$ \\
    {} & $\half + \half\sqrt{3}$ & $\half - \half\sqrt{3}$ \\
  \end{tabular}
\end{table}

\section{The implementation}
A general, implicit Runge-Kutta method (RK method) involves solving a system of non-linear equations every time step.
The exact shape of this system depends on the number of equations in the system of ODEs \emph{and} the number of stages in the method.
If the ODE has $N_{ode}$ equations and the RK method has $N_s$ stages then the total non-linear system has $N_{ode} \times N_s$ equations.
To ease implementation of a general implicit RK method we deduce the general form of the system of equations to solve for a general number of stages and ODEs.
If we denote $\bvec{f}(t,\bvec{y}) = (\bvec{f}_0, \bvec{f}_1, \hdots, \bvec{f}_{N-1})^T$ then we have
\begin{align*}
  \bvec{F} =
  \begin{pmatrix}
    \bvec{f}\left( t + c_0\Delta t, \bvec{y} + \Delta t \sum_{j=0}^{N_s-1} \alpha_{0,j} \bvec{k}_j \right) - \bvec{k}_0 \\
    \vdots \\
    \bvec{f}\left( t + c_0\Delta t, \bvec{y} + \Delta t \sum_{j=0}^{N_s-1} \alpha_{n,j} \bvec{k}_j \right) - \bvec{k}_n \\
    \vdots \\
        \bvec{f}\left( t + c_0\Delta t, \bvec{y} + \Delta t \sum_{j=0}^{N_s-1} \alpha_{N_s-1,j} \bvec{k}_j \right) - \bvec{k}_{N_s-1}
      \end{pmatrix} =
  \bvec{0}
\end{align*}
To solve such a system we employ Newton iteration.
Although Newton iteration is typically costly, we shall see we only need to evaluate $N_{s}$ evaluations of $\bvec{f}$ and its Jacobi matrix $\bvec{J}.$
To find the general expression of the Jacobi matrix, we will apply a general three-stage method to the following system of ODEs:
\begin{align*}
  \bvec{f}(t, \bvec{y}) = -\lambda \bvec{y}
\end{align*}
Then we have
\begin{align*}
  \bvec{F} = \begin{pmatrix}
    -\lambda \bvec{y} - \lambda \Delta t \left( \alpha_{0,0}\bvec{k}_0 + \alpha_{0,1}\bvec{k}_1 + \alpha_{0,2}\bvec{k}_2  \right) - \bvec{k}_0 \\
    -\lambda \bvec{y} - \lambda \Delta t \left( \alpha_{1,0}\bvec{k}_0 + \alpha_{1,1}\bvec{k}_1 + \alpha_{1,2}\bvec{k}_2  \right) - \bvec{k}_1 \\
    -\lambda \bvec{y} - \lambda \Delta t \left( \alpha_{2,0}\bvec{k}_0 + \alpha_{2,1}\bvec{k}_1 + \alpha_{2,2}\bvec{k}_2  \right) - \bvec{k}_2
\end{pmatrix}\end{align*}
The Jacobi matrix of $\bvec{f}$ is given by $-\lambda \bvec{I},$ so the Jacobi matrix of $F$ with respect to $\bvec{k}_i$ becomes
\begin{equation*}
  \bvec{J}_\bvec{k} = \begin{pmatrix}
    -\lambda \Delta t \alpha_{0,0}\bvec{I} - \bvec{I} & -\lambda \Delta t \alpha_{0,1}\bvec{I} & -\lambda \Delta t \alpha_{0,2}\bvec{I} \\
    -\lambda \Delta t \alpha_{1,0}\bvec{I}  & -\lambda \Delta t \alpha_{1,1}\bvec{I}  - \bvec{I} & -\lambda \Delta t \alpha_{1,2}\bvec{I} \\
    -\lambda \Delta t \alpha_{2,0}\bvec{I} & -\lambda \Delta t \alpha_{2,1}\bvec{I} & -\lambda \Delta t \alpha_{2,2}\bvec{I} - \bvec{I} \\
  \end{pmatrix}
\end{equation*}
Now note that $-\lambda \bvec{I}$ is really the Jacobi matrix of $\bvec{J}$ evaluated at specific points. If we write $\bvec{J}\left( t + c_i \Delta t, \bvec{y}_n + \Delta t \sum_{j=0}^{N_s-1} \alpha_{i,j} \bvec{k}_j \right) := \bvec{J}_{i}$ then we have in general that
\begin{align*}
  \bvec{J}_\bvec{k} =& \Delta t\begin{pmatrix}
      \alpha_{0,0} \bvec{J}_0 &  \alpha_{0,1} \bvec{J}_0 &  \alpha_{0,2} \bvec{J}_0 \\
     \alpha_{1,0} \bvec{J}_1  &  \alpha_{1,1} \bvec{J}_1   &  \alpha_{1,2} \bvec{J}_2 \\
     \alpha_{2,0} \bvec{J}_2 &  \alpha_{2,1} \bvec{J}_2 &  \alpha_{2,2} \bvec{J}_2  \\
   \end{pmatrix} - \bvec{I}
\end{align*}
Note that we only need $N_s$ evaluations of the Jacobi matrix of $f$ and not $N_s^2.$
If we have a more general system of ODEs
\begin{equation*}
  \bvec{f}(t,\bvec{y}) = \begin{pmatrix} f_1( t, \bvec{y} ) \\
    f_2( t, \bvec{y} )
  \end{pmatrix}, \qquad \left(\frac{\partial f_k}{\partial z}\right)_i := \frac{ \partial f}{\partial z}\left(t + c_i \Delta t, \bvec{y} + \Delta t \sum_{j=0}^{N_s-1}\alpha_{i,j} \bvec{k}_j\right)
\end{equation*}
with $z = x,y$ and $k = 1,2$ then the total Jacobi matrix system becomes
\begin{align*}
  \bvec{J}_\bvec{k} =& \Delta t\begin{pmatrix}
    \alpha_{0,0} \left( \frac{\partial f_1}{\partial x} \right)_0 & \alpha_{0,0} \left( \frac{\partial f_1}{\partial y} \right)_0 &  \alpha_{0,1} \left( \frac{\partial f_1}{\partial x} \right)_0 & \alpha_{0,1} \left( \frac{\partial f_1}{\partial y} \right)_0 &  \alpha_{0,2} \left( \frac{\partial f_1}{\partial x} \right)_0 & \alpha_{0,2} \left( \frac{\partial f_1}{\partial y} \right)_0 \\
    \alpha_{0,0} \left( \frac{\partial f_2}{\partial x} \right)_0 & \alpha_{0,0} \left( \frac{\partial f_2}{\partial y} \right)_0 &  \alpha_{0,1} \left( \frac{\partial f_2}{\partial x} \right)_0 & \alpha_{0,1} \left( \frac{\partial f_2}{\partial y} \right)_0 &  \alpha_{0,2} \left( \frac{\partial f_2}{\partial x} \right)_0 & \alpha_{0,2} \left( \frac{\partial f_2}{\partial y} \right)_0 \\
    \alpha_{1,0} \left( \frac{\partial f_1}{\partial x} \right)_1 & \alpha_{1,0} \left( \frac{\partial f_1}{\partial y} \right)_1 &  \alpha_{1,1} \left( \frac{\partial f_1}{\partial x} \right)_1 & \alpha_{1,1} \left( \frac{\partial f_1}{\partial y} \right)_1 &  \alpha_{1,2} \left( \frac{\partial f_1}{\partial x} \right)_1 & \alpha_{1,2} \left( \frac{\partial f_1}{\partial y} \right)_1 \\
    \alpha_{1,0} \left( \frac{\partial f_2}{\partial x} \right)_1 & \alpha_{1,0} \left( \frac{\partial f_2}{\partial y} \right)_1 &  \alpha_{1,1} \left( \frac{\partial f_2}{\partial x} \right)_1 & \alpha_{1,1} \left( \frac{\partial f_2}{\partial y} \right)_1 &  \alpha_{1,2} \left( \frac{\partial f_2}{\partial x} \right)_1 & \alpha_{1,2} \left( \frac{\partial f_2}{\partial y} \right)_1 \\
    \alpha_{2,0} \left( \frac{\partial f_1}{\partial x} \right)_2 & \alpha_{2,0} \left( \frac{\partial f_1}{\partial y} \right)_2 &  \alpha_{2,1} \left( \frac{\partial f_1}{\partial x} \right)_2 & \alpha_{2,1} \left( \frac{\partial f_1}{\partial y} \right)_2 &  \alpha_{2,2} \left( \frac{\partial f_1}{\partial x} \right)_2 & \alpha_{2,2} \left( \frac{\partial f_1}{\partial y} \right)_2 \\
    \alpha_{2,0} \left( \frac{\partial f_2}{\partial x} \right)_2 & \alpha_{2,0} \left( \frac{\partial f_2}{\partial y} \right)_2 &  \alpha_{2,1} \left( \frac{\partial f_2}{\partial x} \right)_2 & \alpha_{2,1} \left( \frac{\partial f_2}{\partial y} \right)_2 &  \alpha_{2,2} \left( \frac{\partial f_2}{\partial x} \right)_2 & \alpha_{2,2} \left( \frac{\partial f_2}{\partial y} \right)_2 \\
   \end{pmatrix} - \bvec{I}
\end{align*}
which, in ``prettier'' notation, is
\begin{equation*}
  \bvec{J}_\bvec{k} = \Delta t \begin{pmatrix}
    \alpha_{0,0}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_0 & \alpha_{0,1}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_0 & \alpha_{0,2}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_0 \\
    \alpha_{1,0}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_1 & \alpha_{1,1}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_1 & \alpha_{1,2}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_1 \\
    \alpha_{2,0}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_2 & \alpha_{2,1}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_2 & \alpha_{2,2}\left( \frac{\partial \bvec{f}}{\partial \bvec{y}} \right)_2
  \end{pmatrix} - \bvec{I}
\end{equation*}

Since we have now deduced a general form for the non-linear function that defines all the stages as well as its Jacobi matrix, we can use Newton iteration to solve the system of equations.
For explicit methods, we will of course not resort to Newton iteration because the stages can be computed in a single step.

Instead of calculating the stages using the aforementioned formula, we calculate the following related quantities as per Hairer, N\o rsett and Wanner (Solving Differential Equations II):
\begin{align*}
  \bvec{Y} := \begin{pmatrix}
    \bvec{Y}_1 \\
    \bvec{Y}_2 \\
    \hdots \\
    \bvec{Y}_s
  \end{pmatrix} :=& \Delta t \begin{pmatrix}
    a_{11} \bvec{k}_1 + a_{12}\bvec{k}_2 + \hdots + a_{1s}\bvec{k}_s \\
    a_{21} \bvec{k}_1 + a_{22}\bvec{k}_2 + \hdots + a_{2s}\bvec{k}_s \\
    \vdots \\
    a_{s1} \bvec{k}_1 + a_{s2}\bvec{k}_2 + \hdots + a_{ss}\bvec{k}_s \\
  \end{pmatrix} \\
  =& \Delta t (\bvec{A}\otimes \bvec{I}) ( \bvec{k}_1, \bvec{k}_2, \hdots, \bvec{k}_s)^T \\
  =& \Delta t \begin{pmatrix}
    a_{11}\bvec{I}  & a_{12}\bvec{I} & \hdots & a_{1s}\bvec{I} \\
    a_{21}\bvec{I}  & a_{22}\bvec{I} & \hdots & a_{2s}\bvec{I} \\
    \vdots & \vdots & \ddots & \vdots \\
    a_{s1}\bvec{I}  & a_{s2}\bvec{I} & \hdots & a_{ss}\bvec{I}
  \end{pmatrix} \begin{pmatrix}
    \bvec{k}_1 \\
    \bvec{k}_2 \\
    \vdots \\
    \bvec{k}_s
  \end{pmatrix}
\end{align*}
Since $\bvec{Y} = \Delta t (\bvec{A}\otimes\bvec{I})\bvec{K}$ we have $\bvec{K} = (\Delta t \bvec{A}\otimes\bvec{I})^{-1}\bvec{Y}$ and hence
\begin{equation}
  (\Delta t \bvec{A}\otimes\bvec{I})^{-1}\bvec{Y} = \begin{pmatrix}
    (\bvec{f}(t+c_1\Delta t, \bvec{y} + \bvec{Y}_1)) \\
    (\bvec{f}(t+c_2\Delta t, \bvec{y} + \bvec{Y}_2)) \\
    \vdots \\
    (\bvec{f}(t+c_s\Delta t, \bvec{y} + \bvec{Y}_s))
  \end{pmatrix}
  \label{eqn:system-solve-1}
\end{equation}
Therefore the non-linear system we have to solve each step is given by
\begin{equation}
  \bvec{R} := \bvec{Y} - \Delta t (\bvec{A}\otimes\bvec{I}) \begin{pmatrix}
    (\bvec{f}(t+c_1\Delta t, \bvec{y} + \bvec{Y}_1)) \\
    (\bvec{f}(t+c_2\Delta t, \bvec{y} + \bvec{Y}_2)) \\
    \vdots \\
    (\bvec{f}(t+c_s\Delta t, \bvec{y} + \bvec{Y}_s))
  \end{pmatrix} = 0
  \label{eqn:system-solve-2}
\end{equation}
and its Jacobi matrix is given by
\begin{equation}
  \bvec{J} := \bvec{I}_{Ns\times Ns} - \Delta t (\bvec{A}\otimes\bvec{I}) \begin{pmatrix}
    \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_1) & \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_2) & \hdots & \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_s) \\
    \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_1) & \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_2) & \hdots & \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_s) \\
    \vdots & \vdots & \ddots & \vdots \\
        \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_1) & \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_2) & \hdots & \bvec{J}_\bvec{f}(\bvec{y} + \bvec{Y}_s)
  \end{pmatrix}
  \label{eqn:system-solve-3},
\end{equation}
where we have omitted the time-dependence of the Jacobi matrix $\bvec{J}_\bvec{f}$ of $\bvec{f}.$

For some equations the success of Newton iteration is not impacted much by the Jacobi matrix and for those
cases we can replace $\bvec{J}_\bvec{f}(t + c_i\Delta t, \bvec{y} + \bvec{Y}_i)$ with simply $\bvec{J}_\bvec{f}(t, \bvec{y}).$

\section{A test case}
A simple way to test the integrators is to apply them to an ODE whose solution is known.
We consider
\begin{equation*}
  \frac{\dd}{\dd t} \bvec{y} = \begin{pmatrix}
    -\alpha & -\omega \\
    \omega & -\alpha \end{pmatrix} \bvec{y}
\end{equation*}
The solution can be composed in terms of the eigenvalues $\lambda_\pm$ of the matrix as
\begin{equation*}
  \bvec{y} = \bvec{C}_1e^{\lambda_+t} + \bvec{C}_2e^{\lambda_-t}, \qquad \lambda_\pm = \left[-\alpha \pm i \omega\right].
\end{equation*}
In other words, we have $\bvec{y} = \left[ \bvec{A}_1\cos(\omega t) + \bvec{A}_2 \sin(\omega t) \right]e^{-\alpha t}.$
To find $\bvec{A}_1$ and $\bvec{A}_2$, we apply the initial value $(x,y) = (1,0)^T$ and find
\begin{equation*}
  \bvec{y} = \left[ \begin{pmatrix} 1 \\ 0 \end{pmatrix} \cos(\omega t) +
  \begin{pmatrix} 0 \\ 1 \end{pmatrix} \sin( \omega t ) \right]e^{-\alpha t}
\end{equation*}

\section{Stability regions}

For a given Runge-Kutta method, the stability region can be determined by applying the method to a
simple test problem $y' = \lambda y.$
For a general Runge-Kutta method, the stages are implicitly defined as
\begin{align*}
  \begin{pmatrix} \bvec{k}_0 \\
    \bvec{k}_1 \\
    \vdots \\
    \bvec{k}_{N-1}
  \end{pmatrix}
  =& \lambda \begin{pmatrix} \bvec{y}_0 \\
    \bvec{y}_0 \\
    \vdots \\
    \bvec{y}_0
  \end{pmatrix}
  + \lambda \Delta t \sum_{j=0}^{N-1} \begin{pmatrix}
    \alpha_{0j}\bvec{k}_j \\
    \alpha_{1j}\bvec{k}_j \\
    \vdots \\
    \alpha_{(N-1)j}\bvec{k}_j
  \end{pmatrix} \\
  =& \lambda \begin{pmatrix} \bvec{y}_0 \\
    \bvec{y}_0 \\
    \vdots \\
    \bvec{y}_0
  \end{pmatrix} + \lambda \Delta t \begin{pmatrix}
    \alpha_{00}\bvec{I} &\alpha_{01}\bvec{I} &\hdots &\alpha_{0(N-1)}\bvec{I} \\
    \alpha_{00}\bvec{I} &\alpha_{1j}\bvec{I} &\hdots &\alpha_{1(N-1)}\bvec{I} \\
    \vdots & \vdots & \ddots & \vdots \\
    \alpha_{00}\bvec{I} &\alpha_{(N-1)j}\bvec{I} &\hdots &\alpha_{(N-1)(N-1)}\bvec{I} \\
  \end{pmatrix} \begin{pmatrix}
    \bvec{k}_0 \\
    \bvec{k}_1 \\
    \vdots \\
    \bvec{k}_{N-1}
  \end{pmatrix}
\end{align*}
and hence are given by
\begin{align*}
  \left[ \bvec{I}_{MN\times MN} - \lambda \Delta t \begin{pmatrix}
    \alpha_{00}\bvec{I} &\alpha_{01}\bvec{I} &\hdots &\alpha_{0(N-1)}\bvec{I} \\
    \alpha_{00}\bvec{I} &\alpha_{1j}\bvec{I} &\hdots &\alpha_{1(N-1)}\bvec{I} \\
    \vdots & \vdots & \ddots & \vdots \\
    \alpha_{00}\bvec{I} &\alpha_{(N-1)j}\bvec{I} &\hdots &\alpha_{(N-1)(N-1)}\bvec{I} \\
  \end{pmatrix}\right] \begin{pmatrix} \bvec{k}_0 \\
    \bvec{k}_1 \\
    \vdots \\
    \bvec{k}_{N-1}
  \end{pmatrix} = \lambda
  \begin{pmatrix}
    \bvec{y}_0 \\
    \bvec{y}_0 \\
    \vdots \\
    \bvec{y}_0 \\
  \end{pmatrix},
\end{align*}
where $\bvec{I}_{MN\times MN}$ is a $MN$ by $MN$ identity matrix, as opposed to $\bvec{I}$, which is an $N\times N$ identity matrix, with $M$ the number of equations (1 for our problem) and $N$ the number of stages.
From the solution of this system the new value of $\bvec{y}$ is computed as
\begin{equation*}
  \bvec{y}_1 = \bvec{y}_0 + \Delta t
  \begin{pmatrix}
    \vert & \vert & \hdots & \vert \\
    \bvec{k}_0 & \bvec{k}_1 & \hdots & \bvec{k}_{N-1}\\
        \vert & \vert & \hdots & \vert \\
  \end{pmatrix}\begin{pmatrix}
    b_0 \\
    b_1 \\
    \vdots \\
    b_{N-1}
  \end{pmatrix}
\end{equation*}.
Since we only have one equation, the solution for $\bvec{K} := \left(k_0, k_1, \hdots, k_{N-1}\right)^T$ becomes
\begin{equation*}
  \bvec{K} = \lambda \left[ \bvec{I} - \lambda \Delta t \bvec{A}  \right]^{-1}
  \begin{pmatrix}
    y_0 \\
    y_0 \\
    \vdots \\
    y_0
  \end{pmatrix}
\end{equation*}
and so
\begin{equation*}
  R(\lambda \Delta t) = y_1/y_0 = 1 + \lambda \Delta t \left(b_0, b_1, \hdots, b_{N-1}\right)\left[\bvec{I} - \lambda \Delta t \bvec{A}\right]^{-1}\left(1, 1, \hdots, 1\right)^T
\end{equation*}


\section{Adaptive time step control and embedding methods}

Given a Runge-Kutta method, it is sometimes possible to find a second set of weights $\hat{\bvec{b}}$ in such a way that the approximation $\hat{\bvec{y}} = \bvec{y}_n + \Delta t_n \sum_{i=1}^{N_S} \hat{b}_i\bvec{k}_i$ is also a valid implicit Runge-Kutta method, but with a different order of consistency.
Then, the difference $\bvec{y} - \hat{\bvec{y}}$ is an estimate for the error, since if the original method is of order $m_1$ and the embedded method of order $m_2 < m_1$, then we have
\begin{align*}
  \bvec{y} - \hat{\bvec{y}} =& \bvec{y}_{n+1}^{exact} + \epsilon_1 (\Delta t)^{m_1} - \bvec{y}_{n+1}^{exact} - \epsilon_2 (\Delta t)^{m_2} + \mathcal{O}((\Delta t)^{m_2+1}) \\
  =& \epsilon_{2}(\Delta t)^{m_2} + \mathcal{O}((\Delta t)^{m_2+1})
\end{align*}
in which we recognize that $\bvec{y}-\hat{\bvec{y}}$ is an $m_2$-order consistent approximation to the error in $\hat{\bvec{y}}.$

However, for this type of scheme to work for arbitrary step sizes, it is important that the method given by the coefficients $\hat{\bvec{b}}$ is also L-stable or A-stable. In general, we would like it to have stability characteristics similar to the original method.
To guarantee this we need to go over the same stability checks as the original method.

For example, for the three-stage Lobatto-IIIC method we have
\begin{table}[h]
  \centering
  \begin{tabular}{c|ccc}
    0 & $\frac{1}{6}$ & $-\frac{1}{3}$ & $\frac{1}{6}$ \\
    $\half$ & $\frac{1}{6}$ & $\frac{5}{12}$ & $-\frac{1}{12}$ \\
    1 & $\frac{1}{6}$ & $\frac{2}{3}$ & $\frac{1}{6}$ \\
\hline
    {} & $\frac{1}{6}$ & $\frac{2}{3}$ & $\frac{1}{6}$ \\
  \end{tabular}
\end{table}
from which we can find that
\begin{equation*}
  \bvec{I} - \lambda \Delta t \bvec{A} =
  \begin{pmatrix}
    1 - \lambda \Delta t/6 & \lambda \Delta t/3 & -\lambda \Delta t/6 \\
    -\lambda \Delta t/6 & 1 - 5\lambda \Delta t/12 & \lambda \Delta t/12 \\
    -\lambda \Delta t/6 & -2\lambda \Delta t/3 & 1 - \lambda \Delta t/6
  \end{pmatrix}
\end{equation*}
Now, we are interested in the quantity $\bvec{b}^T(\bvec{I}-z\bvec{A})^{-1}\bvec{e},$ where $z:=\lambda\Delta t$ and $\bvec{e} := (1, 1, 1).$
However, directly computing $(\bvec{I}-z\bvec{A})^{-1}$ is tedious, so we instead determine the $LU$-decomposition first:\\
\begin{align*}
  \bvec{LU} :=&
  \begin{pmatrix}
    1 & 0 & 0 \\
    l_{21} & 1 & 0 \\
    l_{31} & l_{32} & 1
  \end{pmatrix}
  \begin{pmatrix}
    u_{11} & u_{12} & u_{13} \\
    0 & u_{22} & u_{23} \\
    0 & 0 & u_{33}
  \end{pmatrix} \\
  =&\begin{pmatrix}
    u_{11} & u_{12} & u_{13} \\
    l_{21}u_{11} & l_{21}u_{12} + u_{22} & l_{21}u_{13} + u_{23} \\
    l_{31}u_{11} & l_{31}u_{12} + l_{32}u_{22} & l_{31}u_{13} + l_{32}u_{23} + u_{33}
  \end{pmatrix}
\end{align*}
We immediately find the top row of $\bvec{U}$ and first column of $\bvec{L}:$
\begin{align*}
  u_{11} &= 1 - z/6 \qquad u_{12} = z/3 \qquad u_{13} = -z/6 \\
  l_{21} &= -\frac{z/6}{1-z/6} \qquad l_{31} = -\frac{z/6}{1 - z/6}
\end{align*}
To find $u_{22}$ we consider the center equation and find
\begin{equation*}
  -\frac{z^2/18}{1-z/6} + u_{22} = 1 - 5z/12 = \frac{(1 - \frac{5}{12}z)(1 - z/6)}{1 - z/6}
\end{equation*}
from which we isolate
\begin{equation*}
  u_{22} = (1 - z/6)^{-1}\left(1 - \frac{7}{12}z + \frac{5}{72}z^2 + \frac{1}{18}z^2\right) = \frac{1 - \frac{7}{12}z + \frac{1}{8}z^2}{1 - z/6}
\end{equation*}
Now we can substitute this into the third equation of the second column to obtain $l_{32}:$
\begin{align*}
  l_{31}u_{12} + l_{32}u_{22} =& \frac{-z^2/18}{1 - z/6} + l_{32}\frac{1 - \frac{7}{12}z + \frac{1}{8}z^2}{1 - z/6} = -\frac{2z}{3} \\
  l_{32}\frac{1 - \frac{7}{12}z + \frac{1}{8}z^2}{1 - z/6} =& -\frac{2z}{3} + \frac{z^2/18}{1-z/6} = \frac{-2z/3 + 2z^2/18 + z^2/18}{1-z/6} \\
  \rightarrow l_{32} =& \frac{1 - z/6}{1 - \frac{7}{12}z + \frac{1}{8}z^2}\frac{-2z/3 + z^2/6}{1 - z/6} = -\frac{2z}{3}\frac{1 - z/4}{1 - \frac{7}{12}z + \frac{1}{8}z^2}
\end{align*}
To verify correctness, we pause at this point and recompute the first two columns of the product $\bvec{LU}.$ The top rows are trivial, while for the other components we find
\begin{align*}
  l_{21}u_{11} =& -z/6, \qquad l_{31}u_{11} = -z/6\\
  l_{21}u_{12} + u_{22} =& \frac{-z/6}{1-z/6}(z/3) + \frac{1 - \frac{7}{12}z + \frac{1}{8}z^2}{1-z/6} = \frac{1 - \frac{7}{12}z + \frac{1}{8}z^2 - \frac{1}{18}z^2}{1-z/6}
\end{align*}
Since $\frac{1}{8}=\frac{18}{12^2}$ and $\frac{1}{18}=\frac{8}{12^2}$ we end up with
\begin{equation*}
  l_{21}u_{12} + u_{22} = \frac{1 - \frac{7}{12}z + \frac{10}{12^2}z^2}{1-z/6}.
\end{equation*}
The numerator factorizes into $(1-z/6)(1-5z/12)$ so we indeed have $l_{21}u_{12} + u_{22} = 1 - 5z/12.$
Similarly, for the last expression of column 2 we find
\begin{align*}
  l_{31}u_{12} + l_{32}u_{22} =& -\frac{z^2/18}{1-z/6} - \frac{2z}{3}\frac{1 - z/4}{1 - \frac{7}{12}z + \frac{1}{8}z^2}\frac{1 - \frac{7}{12}z + \frac{1}{8}z^2}{1 - z/6} \\
  =& \frac{-z^2/18 - 2z/3 + z^2/6}{1-z/6} = \frac{-2z/3 -z^2/9}{1 - z/6} = -\frac{2z}{3}\frac{1-z/6}{1 - z/6} = -2z/3,
\end{align*}
reassuring us that we made no mistake in any of the terms so far.
The last two terms to find are $u_{23}$ and $u_{33}.$
We have
\begin{align*}
  l_{21}u_{13}  + u_{23} =& \frac{+(z/6)^2}{1-z/6} + u_{23} = \frac{z}{12} \\
  u_{23} =& \frac{z}{12} - \frac{z^2/36}{1-z/6}  = \frac{z/12 - z^2/72 - z^2/36}{1-z/6} = \frac{z}{12}\frac{1-z/2}{1-z/6}
\end{align*}
and finally to find $u_{33}$ we solve
\begin{align*}
  l_{31}u_{13} + l_{32}u_{23} + u_{33} =& \frac{z^2/36}{1 - z/6} - \frac{2z}{3}\frac{z}{12} \frac{1 - z/4}{1 - \frac{7}{12}z + \frac{1}{8}z^2}\frac{1 - z/2}{1 - z/6} + u_{33} = 1 - z/6 \\
  u_{33} =& 1 - \frac{z}{6} - \frac{z^2/36}{1-z/6} + \frac{2z^2}{36}\frac{(1 - z/4)(1-z/2)}{(1-z/6)(1 - \frac{7}{12}z + \frac{1}{8}z^2)}
\end{align*}
The easiest way forward is to combine everything into one fraction. The numerator becomes
\begin{align*}
  & \frac{2z^2}{36}(1-z/4)(1-z/2) + (1 - \frac{z}{6})^2(1-\frac{7}{12}z + \frac{1}{8}z^2) - \frac{z^2}{36}(1-\frac{7}{12}z + \frac{1}{8}z^2) \\
  =& \frac{z^2}{18}\left(1 - \frac{3z}{4} + \frac{z^2}{8}\right) + \left(1 - \frac{z}{3} + \frac{z^2}{36} - \frac{z^2}{36}\right)\left(1 - \frac{7}{12}z + \frac{z^2}{8}\right)  \\
  =& \frac{z^2}{18} - \frac{3z^3}{72} + \frac{z^4}{12^2} + \left(1 - \frac{7}{12}z + \frac{z^2}{8}\right) - \frac{z}{3}\left(1 - \frac{7}{12}z + \frac{z^2}{8}\right) \\
  =& \frac{z^2}{18} - \frac{3z^3}{72} + \frac{z^4}{12^2} + 1 - \frac{7}{12}z + \frac{z^2}{8} - \frac{z}{3} + \frac{7z^2}{36} - \frac{z^3}{24} \\
  =& 1 + z\left(-\frac{7}{12} - \frac{1}{3}\right) + z^2\left(\frac{1}{18} + \frac{1}{8} + \frac{7}{36}\right) + z^3\left(-\frac{3}{72} - \frac{1}{24}\right) + \frac{z^4}{12^2} \\
  =& 1 - \frac{11}{12}z + \frac{3}{8}z^2 - \frac{1}{12}z^3 + \frac{1}{12^2}z^4
\end{align*}
This expression has two real roots, one at $z=6$ and one at around $z\approx 2.6.$
We factorize it as $(1 - z/6)(6z^3 - 36z^2 + 108z - 24)$ so we can cancel one term against the denominator and we finally get
\begin{equation*}
  u_{33} = \frac{6z^3 - 36z^2 + 108z - 24}{1 - \frac{7}{12}z + \frac{1}{8}z^2}
\end{equation*}

Now that we found all coefficients, we can again make sure they are correct. The final equation is
\begin{align*}
  &l_{31}u_{13} + l_{32}u_{23} + u_{33} \\
  &= \frac{+(z/6)^2}{1-z/6} - \frac{2z}{3}\frac{1-z/4}{1 - \frac{7}{12}z + \frac{1}{8}z^2} \frac{z}{12}\frac{1-z/2}{1-z/6} + \frac{6z^3 - 36z^2 + 108z - 24}{1 - \frac{7}{12}z + \frac{1}{8}z^2}
\end{align*}
We combine everything again under the same denominator of $(1-\frac{7}{12}z+\frac{1}{8}z^2)(1-z/6)$


From this we get the full expression for the stability region:
\begin{align*}
  &1 + \frac{z (b_1, b_2, b_3)}{1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24}}\cdot   \begin{pmatrix}
    \frac{z^2}{8} - \frac{7z}{12}+1 & \frac{z^2}{6} - \frac{z}{3} & \frac{z}{6} - \frac{z^2}{24} \\
    \frac{z}{6} - \frac{z^2}{24} & -\frac{z}{3} + 1 & \frac{z^2}{24} - \frac{z}{12} \\
                                                                                           \frac{z^2}{24} + \frac{z}{6} & \frac{2z}{3} - \frac{z^2}{6} & \frac{z^2}{8} - \frac{7z}{12} + 1
                                                                                         \end{pmatrix}
  \begin{pmatrix}
    1 \\
    1 \\
    1
  \end{pmatrix} \\
  =&
  1 + \frac{z (b_1, b_2, b_3)}{1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24}}
  \begin{pmatrix}
    \frac{z^2}{8}-\frac{7z}{12}+1 + \frac{z^2}{6} - \frac{z}{3} + \frac{z}{6} - \frac{z^2}{24} \\
    \frac{z}{6} - \frac{z^2}{24} - \frac{z}{3} + 1 + \frac{z^2}{24} - \frac{z}{12} \\
    \frac{z^2}{24} + \frac{z}{6} + \frac{2z}{3} - \frac{z^2}{6} + \frac{z^2}{8} - \frac{7z}{12} + 1
  \end{pmatrix} \\
  =&
   1 + \frac{z (b_1, b_2, b_3)}{1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24}}
   \begin{pmatrix}
     1 + \frac{3z^2}{24}-\frac{7z}{12} + \frac{4z^2}{24} - \frac{4z}{12} + \frac{2z}{12} - \frac{z^2}{24} \\
     1 + \frac{2z}{12} - \frac{z^2}{24} - \frac{4z}{12} + \frac{z^2}{24} - \frac{z}{12} \\
     1 + \frac{z^2}{24} + \frac{2z}{12} + \frac{8z}{12} - \frac{4z^2}{24} + \frac{3z^2}{24} - \frac{7z}{12}
   \end{pmatrix} \\
     =&
   1 + \frac{z (b_1, b_2, b_3)}{1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24}}
   \begin{pmatrix}
     1 - \frac{9}{12} z + \frac{6}{24}z^2 \\
     1 - \frac{3}{12} z \\
     1 - \frac{1}{12} z
   \end{pmatrix}
\end{align*}
Further simplification leads to
\begin{align*}
  1 + \frac{z (b_1 + b_2 + b_3) + \frac{z^2}{12}\left(-9b_1 - 3b_2 - b_3\right) +  b_1\frac{6z^3}{24}}{1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24}}
\end{align*}
For the given method we have $b_1 = b_3 = 1/6$ and $b_2 = 2/3 = 4/6$ so we obtain
\begin{equation*}
  \left(1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24}\right)^{-1}\left(1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24} + z - \frac{22 z^2}{12} + \frac{z^3}{24} \right)
\end{equation*}
which finally reduces to
\begin{equation*}
  \frac{1 + \frac{z}{4} - \frac{1}{12}z^2}{1 - \frac{3z}{4} + \frac{z^2}{4} - \frac{z^3}{24}}
\end{equation*}
From this expression, we see that the method is L-stable, because in the limit of $z\downarrow -\infty$, clearly the denominator dominates and we obtain 0.
To determine A-stability, we need to find the contours where the stability function is equal to 1.

Even this simple 3-stage method was a pain, so it would be great if there was a more straightforward way to go about determining the stability region of the method. In particular, it would be helpful to be able to find $\bvec{b}$ so that the method is L-stable or at least A-stable.

If the absolute value of all eigenvalues of $\bvec{A}$ are less than 1, we can for example use the identity
\begin{equation*}
  (\bvec{I}-z\bvec{A})^{-1} = \sum_{n=0}^{\infty} z^n \bvec{A}^n
\end{equation*}
If we decompose $\bvec{A}$ into its eigenvalues and eigenvectors this furthermore becomes
\begin{equation*}
  (\bvec{I}-z\bvec{A})^{-1} = \sum_{n=0}^{\infty} z^n \bvec{P}\bvec{D}^n \bvec{P}^{-1}
\end{equation*}
where $\bvec{P}$ is a matrix containing the eigenvectors of $\bvec{A}$ and $\bvec{D}$ is a diagonal matrix with the eigenvalues of $\bvec{A}$ on its diagonal.

Multiplying left and right with $\bvec{b}$ and $\bvec{e}$ gives us
\begin{equation*}
  \bvec{b}^T(\bvec{I}-z\bvec{A})^{-1}\bvec{e} = \bvec{b}^T\bvec{P}  \sum_{n=0}^{\infty} z^n \bvec{D}^n  (\bvec{P}^{-1}\bvec{e})
\end{equation*}
This means we can in principle just evaluate $\bvec{b}^{T}\bvec{P}$ once and $\bvec{P}^{-1}\bvec{e}$ and we can generate the full stability polynomial.

As proof of concept, we consider the two-stage Lobatto-IIIC method with
\begin{equation*}
  \bvec{A} = \begin{pmatrix}
    \half & -\half \\
    \half & \half
  \end{pmatrix}
\end{equation*}
This matrix has characteristic polynomial $\rho(\lambda) = (\lambda - \half)^2 + \frac{1}{4} = \lambda^2 - \lambda + \half.$ from which we find eigenvalues $\lambda = (1 \pm i)/2.$ To verify, $(1 \pm i)^2/4 = (1 \pm 2i + i^2)/4 = \pm i/2.$ Therefore, substituting the eigenvalue back in $\rho(\lambda)$ gives us $\pm i/2 - \half \mp i/2 + \half=0$ so we indeed have the proper roots.

To find the eigenvectors, we solve
\begin{equation*}
  \begin{pmatrix}
    \frac{1 \pm i}{2} - \half & \half \\
    -\half & \frac{1 \pm i}{2} - \half
  \end{pmatrix}
  \begin{pmatrix}
    x \\ y
  \end{pmatrix} =
  \begin{pmatrix}
    \pm \frac{i}{2} & \half \\
    -\half & \pm \frac{i}{2}
  \end{pmatrix}
  \begin{pmatrix}
    x \\ y
  \end{pmatrix} =
  \half
  \begin{pmatrix}
    \pm i x + y \\
    -x \pm i y
  \end{pmatrix}
  = \bvec{0}.
\end{equation*}
For the $+$ we find $\frac{1}{\sqrt{2}}(i, 1)^T$ as eigenvector and for $-$ we obtain $\frac{1}{\sqrt{2}}(-i, 1)^T.$ The eigenvectors are already normalized and we find that
\begin{align*}
  \bvec{P}\bvec{P}^{-1} =
  \half \begin{pmatrix}
    i & -i \\
    1 & 1
  \end{pmatrix}
  \begin{pmatrix}
   -i & 1 \\
    i & 1
  \end{pmatrix} = \half
  \begin{pmatrix}
    -i^2 - i^2 & i - i \\
     -i + i  & 1 + 1
  \end{pmatrix} = \bvec{I}
\end{align*}
Therefore, we also know that
\begin{align*}
  \bvec{P}\bvec{D}\bvec{P}^{-1} &=
  \frac{1}{4} \begin{pmatrix}
    i & -i \\
    1 &  1
  \end{pmatrix}
  \begin{pmatrix}
    1 + i & \\
    & 1 - i
  \end{pmatrix}
  \begin{pmatrix}
    -i & 1 \\
    i &  1
  \end{pmatrix} = \frac{1}{4}
  \begin{pmatrix}
    i & -i \\
    1 & 1
  \end{pmatrix}
  \begin{pmatrix}
    -i - i^2 & 1 + i \\
    i - i^2 & 1 - i
  \end{pmatrix} \\
  =& \frac{1}{4}
  \begin{pmatrix}
    i & -i \\
    1 &  1
  \end{pmatrix}
  \begin{pmatrix}
   1 - i & 1 + i \\
   1 + i & 1 - i \\
  \end{pmatrix} = \frac{1}{4}
  \begin{pmatrix}
    i - i^2 - i - i^2 & i + i^2 -i + i^2 \\
    1 - i + 1 + i & 1 + i + 1 - i
  \end{pmatrix} \\
  =& \frac{1}{4}\begin{pmatrix}
    2 & -2 \\
    2 & 2
  \end{pmatrix} =
  \begin{pmatrix}
    \half & -\half \\
    \half & \half
  \end{pmatrix} = \bvec{A}
\end{align*}
Using this decomposition it is also straightforward to determine $(\bvec{I}-z\bvec{A})^{-1},$ since we can write
\begin{align*}
  (\bvec{I} - z\bvec{A})^{-1} =&   (\bvec{P}\bvec{P}^{-1} - z\bvec{P}\bvec{D}\bvec{P}^{-1})^{-1} = \left(\bvec{P}\left(\bvec{I} - z\bvec{D}\right)\bvec{P}^{-1}\right)^{-1} = \\
  =& \bvec{P}\left(\bvec{I} - z\bvec{D}\right)^{-1}\bvec{P}^{-1}
\end{align*}
We already know what $\bvec{D}$ is so we can simply write
\begin{align*}
  (\bvec{I} - z\bvec{D})^{-1}
  \begin{pmatrix}
    1 - \half z(1 + i) & \\
    & 1 - \half z(1 - i)
  \end{pmatrix}^{-1} =
  2\begin{pmatrix}
    \frac{1}{2 - z(1 + i)} & \\
    & \frac{1}{2 - z(1 - i)}
  \end{pmatrix}
\end{align*}
Reconstructing the inverse matrix is now simply evaluating
\begin{align*}
  (\bvec{I} - z\bvec{A})^{-1} =&
  \half\cdot2 \begin{pmatrix}
    i & -i \\
    1 &  1
  \end{pmatrix}
  \begin{pmatrix}
    \frac{1}{2 - z(1 + i)} & \\
    & \frac{1}{2 - z(1 - i)}
  \end{pmatrix}
  \begin{pmatrix}
   -i & 1 \\
    i & 1
  \end{pmatrix}
\end{align*}
but it's even more convenient to directly use the expression for the stability region and compute directly:
\begin{align*}
  R(Z) = y(t_{n+1})/y(t_n) =& 1 + z \bvec{b}\cdot(\bvec{I}-z\bvec{A})^{-1} \bvec{e} \\
  =& 1 + z \bvec{b}\cdot \bvec{P}(\bvec{I}-z\bvec{D})^{-1} \bvec{P^{-1}e}
\end{align*}
We can easily compute the outermost terms to obtain
\begin{align*}
  R(z) =& 1 + z\begin{pmatrix}
    b_1i + b_2 \\
   -b_1i + b_2
  \end{pmatrix} \begin{pmatrix}
    \frac{1}{2 - z(1+i)} & \\
    & \frac{1}{2 - z(1-i)}
  \end{pmatrix}
  \begin{pmatrix}
    1 - i \\
    1 + i
  \end{pmatrix}
\end{align*}
The right-most product becomes $\left((1-i)/(2-z - iz), (1+i)/(2 - z + iz)\right)^T$ and therefore
\begin{align*}
  R(z) =& 1 + z \begin{pmatrix}
    b_2 + b_1i  \\
   b_2 - b_1i
  \end{pmatrix}^T\begin{pmatrix}
    \frac{1-i}{2-z - iz} \\
    \frac{i+1}{2 - z + iz}
  \end{pmatrix} = 1 + z \left[ \frac{(b_2+ib_1)(1-i)}{2 -z - iz} + \frac{(b_2- ib_1)(1+i)}{2-z+iz}\right] \\
  =& 1 + \frac{z[(b_2 + ib_1 - ib_2 + b_1)(2 - z + iz) + (b_2 - ib_1 + ib_2 + b_1)(2 - z - iz)]}{(2-z)^2 + z^2} \\
  =& 1 + \frac{z[(b_1 + b_2 + i(b_1 - b_2))(2 - z + iz) + (b_2(1 + i) + b_1(1 - i))(2 - z - iz)]}{(2-z)^2 + z^2} \\
  =& 1 + \frac{z(2-z)\left[b_1 + b_2 + ib_1 - ib_2 + b_2 + ib_2 + b_1 - ib_1\right]}{(2-z)^2 + z^2}  \\
        &+ \frac{iz^2 \left[  (b_1 + b_2 + ib_1 - ib_2 - b_2 - ib_2 - b_1 + ib_1)\right]  }{(2-z)^2 + z^2} \\
  =& 1 + \frac{z(2-z)\left(2b_1 + 2b_2\right) + iz^2\left(2i(b_1 - b_2)\right)}{4 - 4z + 2z^2} \\
  =& 1 + \frac{4(b_1+b_2)z - 2(b_1+b_2)z^2 - 2z^2 (b_1-b_2)}{4 - 4z + 2z^2} = 1 + \frac{4(b_1+b_2)z - 4 b_1 z^2}{4 - 4z + 2z^2} \\
  =& 1 + \frac{(b_1+b_2)z - b_1z^2}{1 - z + \half z^2} = \frac{1 + (b_1+b_2 - 1)z + (\half - b_1) z^2}{1 - z + \half z^2}
\end{align*}
Choosing $b_1 = b_2 = \half$ clearly recovers the Pad\'{e} approximant $P(n-2, n)$ and the method is L-stable.

Applying the method directly leads to the following:
\begin{align*}
  \begin{pmatrix}
    k_1 \\ k_2
  \end{pmatrix} =  \begin{pmatrix} (\lambda y_0 + z (\half k_1 - \half k_2) \\
    (\lambda y_0 + z(\half k_1 + \half k_2))
    \end{pmatrix}
\end{align*}
From the first equation we have $k_1(1 - \half z) = \lambda y_0 - \half z k_2$ and so
\begin{align*}
  k_2 \left(1 - \half z\right)^2 =& \lambda y_0\left(1 - \half z\right) + \half z \left(\lambda y_0 - \half z k_2\right) \\
  k_2\left[(1 - \half z)^2 + \frac{1}{4}z^2\right] =& \lambda y_0 \\
  k_2 =& \frac{\lambda y_0}{1 - z + \half z^2}, \qquad k_1 = \frac{\lambda y_0(1 - z)}{1 - z + \half z^2} = (1-z)k_2
\end{align*}
This means that for the next value for $y$ we have
\begin{align*}
  y_1 = y_0 + \Delta t (b_1k_1 + b_2k_2) = y_0 + z y_0 \frac{b_1(1-z)+b_2}{1 - z + \half z^2}
\end{align*}
and the stability function becomes
\begin{align*}
  R(z) =& y_1/y_0 = 1 + z\frac{b_1(1-z) + b_2}{1 - z + \half z^2} = \frac{1 - z + \half z^2 + (b_1+b_2)z - b_1z^2}{1 - z + \half z^2}
\end{align*}
At $b_1=b_2=\half$ we find
\begin{align*}
  R(z) =&  \frac{1 - z + \half z^2 + z - \half z^2}{1 - z + \half z^2} = \frac{1}{1 - z + \half z^2}
\end{align*}
which is correct since Lobatto IIIC methods of order $n$ correspond to the Pad\'{e} approximant $P(n-2, n).$


We now use the above method involving eigenvalues to derive two sets of coefficients for the 3-stage Lobatto IIIC method. We have the following Butcher tableau:

\begin{table}[h]
  \centering
  \begin{tabular}{c|cccc}
    0 & $\frac{1}{6}$ & $-\frac{1}{3} $ & $\frac{1}{6}$ \\
    $\half$ & $\frac{1}{6}$ & $\frac{5}{12}$ & $-\frac{1}{12}$ \\
    1 & $\frac{1}{6}$ & $\frac{2}{3}$ &  $\frac{1}{6}$ \\ \hline
    {} & $\frac{1}{6}$ & $\frac{2}{3}$ & $\frac{2}{6}$
  \end{tabular}
\end{table}
The inverse is given by
\begin{align*}
  \bvec{A}^{-1} =& \begin{pmatrix}
     3 &  4 & -1 \\
    -1 &  0 &  1 \\
     1 & -4 &  3
  \end{pmatrix} \\
  \bvec{A}^{-1}\bvec{A} =& \frac{1}{12} \begin{pmatrix}
     3 &  4 & -1 \\
    -1 &  0 &  1 \\
     1 & -4 &  3
  \end{pmatrix} \begin{pmatrix}
    2 & -4 & 2 \\
    2 & 5 & -1 \\
    2 & 8 & 2
  \end{pmatrix} \\
  =& \frac{1}{12}\begin{pmatrix}
    6 + 8 - 2 & -12 + 20 - 8 & 6 - 4 - 2 \\
    -2 + 2 & 4 + 8 & -2 + 2 \\
    2 - 8 + 6 & -4 - 20 + 24 & 2 + 4 + 6
  \end{pmatrix} = \begin{pmatrix}
    12 & 0 & 0 \\
    0 & 12 & 0 \\
    0 & 0 & 12
  \end{pmatrix} = \bvec{I}
\end{align*}
We can find the eigenvalues from $\det (\lambda \bvec{I} - \bvec{A})) = 0:$
\begin{align*}
  \det (\lambda \bvec{I} - \bvec{A}) =& \left| \begin{matrix}
    \lambda - \frac{1}{6} & \frac{1}{3} & -\frac{1}{6} \\
    -\frac{1}{6} & \lambda - \frac{5}{12} & \frac{1}{12} \\
    -\frac{1}{6} & -\frac{2}{3} & \lambda - \frac{1}{6}
  \end{matrix}
  \right| = (\lambda - \frac{1}{6}) \left[ (\lambda - \frac{5}{12})(\lambda - \frac{1}{6}) + \frac{1}{18} \right] \\
                                      & + \frac{1}{6}\left[\frac{\lambda}{3} - \frac{1}{18} - \frac{1}{9}\right] - \frac{1}{6}\left[ \frac{1}{36} + \frac{1}{6}(\lambda - \frac{5}{12}) \right] \\
  =& (\lambda - \frac{1}{6})^2(\lambda - \frac{5}{12}) + \frac{1}{18}(\lambda - \frac{1}{6}) + \frac{\lambda}{18} - \frac{1}{36} - \frac{1}{6^3} - \frac{1}{6^2}(\lambda - \frac{5}{12}) \\
  =& (\lambda^2 - \frac{1}{3}\lambda + \frac{1}{36})(\lambda - \frac{5}{12}) + \frac{1}{36}\lambda - \frac{1}{18\cdot 6} - \frac{1}{6^3} + \frac{5}{2\cdot 6^3} \\
  =& \lambda^3 - \frac{1}{3}\lambda^2 + \frac{1}{36}\lambda - \frac{5}{12}\lambda^2 + \frac{5}{36}\lambda - \frac{5}{12\cdot 36} + \frac{1}{36}\lambda - \frac{1}{2\cdot 6^3} \\
  =& \lambda^3 - \frac{3}{4}\lambda^2 + \frac{7}{36}\lambda - \frac{1}{2\cdot 6^3} = 0
\end{align*}
The solution is quite complicated, we use Maxima to obtain the roots:
\begin{align*}
  \lambda =& -ac_+ + \frac{c_-}{48a} + \frac{1}{4}, -ac_- + \frac{c_+}{48a}  + \frac{1}{4}, a - \frac{1}{4a} + \frac{1}{4}, \\
  a =& \left(\frac{1}{96^{3/2}} + \frac{1}{192}\right)^{1/3}, \qquad
c_\pm = \half \pm \sqrt{3}i\half
\end{align*}
This is also not super handy. More convenient is to let Maxima do some work on the direct equation
\begin{equation*}
  R(z) = 1 + z \left(b_0, b_1, b_2\right)\left[\bvec{I} - z \bvec{A}\right]^{-1}\left(1, 1, 1\right)^T
\end{equation*}
We can simplify and expand $[\bvec{I}-z\bvec{A}]^{-1}(1,1,1)^T$ to
\begin{align*}
  [\bvec{I}-z\bvec{A}]^{-1}(1,1,1)^T = \frac{1}{1 - \frac{3}{4}z + \frac{1}{4}z^2 - \frac{z^3}{24}}\begin{pmatrix}
    1 - \frac{3}{4}z + \frac{1}{4}z^2 \\
    1 - \frac{3}{4}z \\
    1 + \frac{3}{4}z
    \end{pmatrix}
\end{align*}
and obtain
\begin{align*}
  &1 + z \bvec{b}\cdot[\bvec{I}-z\bvec{A}]^{-1}(1,1,1)^T \\
  &= \frac{1 - \frac{3}{4}z + \frac{1}{4}z^2 - \frac{z^3}{24} + b_1z(1 - \frac{3}{4}z + \frac{1}{4}z^2) + b_2z(1 - \frac{1}{4}z) + b_3z(1 + \frac{1}{4}z)}{1 - \frac{3}{4}z + \frac{1}{4}z^2 - \frac{z^3}{24}}
\end{align*}
We can collect all terms of equal powers of $z$ in the numerator:
\begin{align*}
  \frac{1 + \left(b_1 + b_2 + b_3 - \frac{3}{4}\right)z + \left(\frac{1}{4} - \frac{3}{4}b_1 - \frac{1}{4}b_2 + \frac{1}{4}b_3 \right)z^2 + \left(\frac{1}{4}b_1 - \frac{1}{24}\right)z^3}{1 - \frac{3}{4}z + \frac{1}{4}z^2 - \frac{1}{24}z^3}
\end{align*}
To recover the Pad\'{e} approximant $P(1, 3)$ we require
\begin{align*}
  b_1+b_2+b_3 - \frac{3}{4} = \frac{1}{4}, \quad 1 + b_3 - b_2 - 3b_1 = 0, \quad b_1 - \frac{1}{6} = 0
\end{align*}
Clearly we obtain $b_1 = 1/6,$ $b_2+b_3 = 5/6$ and $b_3-b_2 = -\frac{1}{2}$ and adding the latter two leads to $2b_3 = \frac{5}{6} - \frac{3}{6} = \frac{1}{3}$ so we find $b_3 = \frac{1}{6}$ and $b_2=\frac{2}{3}.$

For an embedded pair, we would still require $b_1=1/6$ because the numerator needs to be a power lower than $z^3.$ A second requirement is always that $b_1=b_2=b_3=1$ or else the method is not consistent.
This leaves us with only a single unknown, since we have $b_1 = \frac{1}{6}$ and $b_2 = \frac{5}{6} - b_3.$ The 3-stage Lobatto IIIC methods has no $z^2$ term in the numerator, but we can choose an alternate value $b_2$ to keep it.
For example, we can choose $b_2=\half$ and obtain $b_3 = \frac{5}{6} - \frac{3}{6} = \frac{1}{3}.$
This means that for the stability region we have
\begin{equation*}
  R(z) = \frac{1 + \frac{1}{4}z + \frac{1}{12}z^2}{1 - \frac{3}{4}z + \frac{1}{4}z^2 - \frac{1}{24}z^3}
\end{equation*}
Clearly, in the limit of $z\rightarrow -\infty,$ we obtain $R(z) \rightarrow 0.$

We might as well make the choice for $b_2$ optimal in some sense. For example, we could maximize the accuracy/order of consistency of the method, if possible.
The Taylor series expansion of the stability region is
\begin{align*}
  R(z) =& 1 + (b_1+b_2+b_3)z + \half \left(2b_2 + b_3\right)z^2 + \frac{1}{6}\left(\frac{3}{4}b_2 + 6b_3\right)z^3 \\
  &+ \frac{1}{24}\left(b_1 + \frac{1}{4}b_2 + 4b_3\right)z^4 + \mathcal{O}(z^5)
\end{align*}
Substitution of $b_1=\frac{1}{6}$, which is required for L-stability, and $b_2 = \frac{5}{6}-b_3,$ which follows from $b_1+b_2+b_3=1,$ we have
\begin{align*}
  R(z) =& 1 + z + \half \left(\frac{5}{3} - b_3\right)z^2 + \frac{1}{6}\left(\frac{3}{4}\left(\frac{5}{6} - b_3\right) + 6b_3\right)z^3 \\
  &+ \frac{1}{24}\left(\frac{1}{6} + \frac{1}{4}\left(\frac{5}{6} - b_3\right) + 4b_3\right)z^4 + \mathcal{O}(z^5)
\end{align*}
The Taylor series for $\exp(z)$ up to fifth order is of course $1 + z + \half z^2 + \frac{1}{6}z^3 + \frac{1}{24}z^4 + \mathcal{O}(z^5)$ so we require $b_3 = 2/3.$
We summarize the third order Lobatto IIIC method with embedded method in the following Butcher tableau:
\begin{table}[h]
  \centering
  \begin{tabular}{c|ccc}
    $0$ & $\frac{1}{6}$ & $-\frac{1}{3}$ & $\frac{1}{6}$ \\
    $\half$ & $\frac{1}{6}$ & $\frac{5}{12}$ & $-\frac{1}{12}$ \\
    $1$ & $\frac{1}{6}$ & $\frac{2}{3}$  & $\frac{1}{6}$ \\ \hline
    {} & $\frac{1}{6}$ & $\frac{2}{3}$  & $\frac{1}{6}$ \\
    {} & $\frac{1}{6}$ & $\frac{1}{6}  $& $\frac{2}{3}$
  \end{tabular}
\end{table}
Interestingly, the coefficients for the embedded $\bvec{b}$ are a cyclic rotation of the original $b.$



\subsection{A note on stability of embedded methods}
In the above section, we looked at Lobatto IIIC methods and saw that because the main Lobatto IIIC method corresponds to a Pad\'{e} $P(n-2, n)$ approximation, we can choose parameters $\bvec{b}$ for the embedded method such that it is of a lower order but still $L$-stable.
I think that for Radau methods that is, in general, not possible.
Consider the 2-stage Radau IIA method, given by
\begin{table}[h]
  \centering
  \begin{tabular}{c|ccc}
    $\frac{1}{3}$ & $\frac{5}{12}$ & $-\frac{1}{12} $ \\
    $1$ & $\frac{3}{4}$ & $\frac{1}{4}$ \\ \hline
    {} & $\frac{3}{4}$ & $\frac{1}{4}$
  \end{tabular}
\end{table}
so we have
\begin{align*}
  \bvec{I}-z\bvec{A} =& \begin{pmatrix}
    1 - \frac{5}{12}z & \frac{1}{12}z \\
    -\frac{3}{4}z & 1 - \frac{1}{4}z
  \end{pmatrix}, \\
  (\bvec{I}-z\bvec{A})^{-1} =& \frac{1}{(1-\frac{5}{12}z)(1-\frac{1}{4}z) + \frac{1}{16}z^2}
  \begin{pmatrix}
    1 - \frac{1}{4}z & -\frac{1}{12}z \\
    \frac{3}{4}z & 1 - \frac{5}{12}z
  \end{pmatrix}
\end{align*}
% test
% \begin{align*}
%   (\bvec{I}-z\bvec{A})^{-1}(\bvec{I}-z\bvec{A}) =&
%   \frac{1}{(1-\frac{5}{12}z)(1-\frac{1}{4}z) + \frac{1}{16}}
%   \begin{pmatrix}
%     1 - \frac{5}{12}z & -\frac{3}{4} \\
%     \frac{1}{12} & 1 - \frac{1}{4}z
%   \end{pmatrix}
%   \begin{pmatrix}
%     1 - \frac{1}{4}z & \frac{3}{4} \\
%     -\frac{1}{12} & 1 - \frac{5}{12}z
%   \end{pmatrix} \\
%   =& \frac{1}{(1-\frac{5}{12}z)(1-\frac{1}{4}z) + \frac{1}{16}}
%      \begin{pmatrix}
%        (1-\frac{5}{12}z)(1 - \frac{1}{4}z) + \frac{1}{16} & 0 \\
%        0 & \frac{1}{16} + (1-\frac{1}{4}z)(1-\frac{5}{12}z)
%      \end{pmatrix} = \bvec{I}
% \end{align*}

and the stability region is
\begin{align*}
  1 + z (b_0, b_1) (\bvec{I}-z\bvec{A})^{-1}(1,1)^T =& 1 + \frac{z}{(1-\frac{5}{12}z)(1-\frac{1}{4}z) + \frac{1}{16}z^2} (b_0, b_1) \begin{pmatrix}
    1 - \frac{1}{3}z \\
    1 + \frac{1}{3}z
  \end{pmatrix} \\
  =& 1 + z\frac{(b_0 + b_1) - \frac{b_0-b_1}{3}z}{1 - \frac{2}{3}z + \frac{1}{6}z^2}
\end{align*}
Simplifying this to one fraction leads to
\begin{align*}
  R(z) = \frac{1 + \left(b_0 + b_1 - \frac{2}{3}\right)z + \left(\frac{1}{6}-\frac{1}{3}b_0+\frac{1}{3}b_1\right)z^2}{1 - \frac{2}{3}z + \frac{1}{6}z^2}
\end{align*}
Consistency requires $b_0+b_1=1$ and for the Radau IIA method we have $b_1=\frac{1}{4}$ to obtain
\begin{align*}
  R(z) = \frac{1 + \frac{1}{3}z}{1 - \frac{2}{3}z + \frac{1}{6}z^2}
\end{align*}
We see that in this case, because $L$-stability is guaranteed by virtue of $b_0=3b_1=\frac{3}{4}$, we cannot make another choice for $\bvec{b}$ that will guarantee an $L$-stable embedded method that is still consistent. The best we can do is $A$-stability, for which we might not be able to guarantee in general.
The natural choice based on the Pad\'{e} approximation is to have the numerator match the Taylor series of $\exp(\frac{1}{3}z)$ up to second degree:
\begin{equation*}
  1 + \frac{1}{3}z + \left(\frac{1}{6}-\frac{1}{3}b_0+\frac{1}{3}b_1\right)z^2 := 1 + \frac{1}{3}z + \frac{1}{18}z^2
\end{equation*}

which leads to $(b_1 - b_0) = -\frac{1}{3}.$ Combining with $b_0+b_1=1$ we find $2b_1 = \frac{2}{3}$ so $b_1=\frac{1}{3}$ and $b_0 = \frac{2}{3}.$ We end up with
\begin{equation*}
  R(z) = \frac{1 + \frac{1}{3}z + \frac{1}{18}z^2}{1 - \frac{2}{3}z + \frac{1}{6}z^2}
\end{equation*}
In the limit of $z \downarrow -\infty$ we see that $R(z) \rightarrow \frac{1}{3}.$
Picking $b_0=b_1$ is about the worst we can do, it makes the method marginally $A$-stable with $R(z\downarrow -\infty)=1.$
For optimal consistency, we can take the product of the Taylor series to find
\begin{align*}
  &\left[1 + \frac{1}{3}z + \left(\frac{1}{6} - \frac{1}{3}b_0 + \frac{1}{3}b_1\right)z^2\right]\left[1 + \frac{2}{3}z + \frac{1}{6}z^2\right] \\
  =& 1 + z + \half z^2\left(\frac{1}{3} - \frac{2b_0}{3} + \frac{2b_1}{3} + \frac{1}{3} + \frac{4}{9}\right) + \frac{1}{6}z^3 \left(\frac{1}{3}  + \frac{2}{3}\left(1 - 2b_0 + 2b_1\right)\right)  + \mathcal{O}(z^4)
\end{align*}
Since we only have $b_1-b_0$ appearing in this form, we only have one degree of freedom, and choosing second-order consistency leads to
\begin{equation*}
  \frac{10}{9} + \frac{6}{9}(b_1-b_0) = 1
\end{equation*}
from which we obtain $b_1-b_0 = -\frac{1}{6}$ which is consistent with our previous choice of $b_0 = \frac{2}{3} = 2b_1.$
We summarize the third order Radau IIA method with embedded method in the following Butcher tableau:
\begin{table}[h]
  \centering
  \begin{tabular}{c|ccc}
    $\frac{1}{3}$ & $\frac{5}{12}$ & $-\frac{1}{12} $ \\
    $1$ & $\frac{3}{4}$ & $\frac{1}{4}$ \\ \hline
    {} & $\frac{3}{4}$ & $\frac{1}{4}$ \\
    {} & $\frac{1}{3}$ & $\frac{2}{3}$
  \end{tabular}
\end{table}

\subsection{More coefficients}
In the appendix we derive/show more coefficients for embedded IRK methods that Rehuel supports.

\clearpage{}

\section{Multistep methods}
Another approach to solving ODEs is to use multistep methods.
The most sensible methods for solving stiff equations are backward-differentiation formula methods.
The idea is straight-forward: We can integrate an ODE with RHS $\bvec{y}' = f(t, \bvec{y})$ by applying a finite difference method to the LHS. The most straightforward
\begin{equation*}
  \bvec{y}_{n+1} - \bvec{y}_n = \Delta t\cdot f(t_n + \Delta t, \bvec{y}_{n+1}).
\end{equation*}
This is clearly just a first-order finite difference approximation to $\Delta t \cdot f(t_n + \Delta t, \bvec{y}_{n+1}.$
The idea behind BDFs is to replace the right-hand-side with an interpolation polynomial rather than only the time derivative at $n+1$.
We can e.g. use Lagrange polynomials for the last four points. If the points are also equidistant, we have $t_n - t_m=(n-m)\Delta t:$
\begin{align*}
  p_1(t) =& \frac{t - t_n}{t_{n+1}-t_{n}} \frac{t - t_{n-1}}{t_{n+1} - t_{n-1}} \frac{t - t_{n-2}}{t_{n+1} - t_{n-2}} = \frac{(t-t_n)(t-t_{n-1})(t-t_{n-2})}{6(\Delta t)^3} \\
  p_2(t) =& \frac{t - t_{n+1}}{t_{n}-t_{n+1}} \frac{t - t_{n-1}}{t_{n} - t_{n-1}} \frac{t - t_{n-2}}{t_{n} - t_{n-2}} = -\frac{(t - t_{n-1})(t-t_{n-2})(t-t_{n+1})}{2(\Delta t)^3} \\
  p_3(t) =& \frac{t - t_{n+1}}{t_{n-1}-t_{n+1}} \frac{t - t_{n-2}}{t_{n-1} - t_{n-2}} \frac{t - t_{n}}{t_{n-1} - t_{n}} = \frac{(t - t_{n+1})(t-t_{n-2})(t-t_{n})}{2 (\Delta t)^3} \\
  p_4(t) =& -\frac{t - t_{n-1}}{t_{n-2}-t_{n-1}} \frac{t - t_{n}}{t_{n-2} - t_{n}} \frac{t - t_{n+1}}{t_{n-2} - t_{n+1}} = \frac{(t - t_{n-1})(t - t_{n})(t-t_{n+1})}{6 (\Delta t)^3}
\end{align*}
We construct a linear combination of these polynomials to go through $\bvec{y}_{n+1},~\bvec{y}_n,~\bvec{y}_{n-1}$ and $\bvec{y}_{n-2}:$
\begin{equation*}
  P(t) = \bvec{y}_{n+1}p_1(t) + \bvec{y}_{n}p_2(t) + \bvec{y}_{n-1}p_3(t) + \bvec{y}_{n-2}p_4(t).
\end{equation*}
The derivative of this polynomial is
\begin{equation*}
P(t) = \bvec{y}_{n+1}p_1'(t) + \bvec{y}_{n}p_2'(t) + \bvec{y}_{n-1}p_3'(t) + \bvec{y}_{n-2}p_4'(t).
\end{equation*}
The derivatives of each polynomial separately are
\begin{align*}
  p_1'(t) =& \frac{(t-t_{n-1})(t-t_{n-2}) + (t-t_n)(t-t_{n-2}) + (t-t_n)(t-t_{n-1})}{6(\Delta t)^3} \\
  p_2'(t) =& -\frac{(t-t_{n-2})(t-t_{n+1}) + (t - t_{n-1})(t-t_{n+1}) + (t - t_{n-1})(t-t_{n-2})}{2(\Delta t)^3} \\
  p_3'(t) =& \frac{(t-t_{n-2})(t-t_{n}) + (t - t_{n+1})(t-t_{n}) + (t - t_{n+1})(t-t_{n-2})}{2 (\Delta t)^3} \\
  p_4'(t) =& -\frac{(t - t_{n})(t-t_{n+1}) + (t - t_{n-1})(t-t_{n+1}) + (t - t_{n-1})(t - t_{n})}{6 (\Delta t)^3}
\end{align*}
Evaluating them at $t_{n+1}$ allows us to evaluate $P'(t)$ at $t_{n+1}:$
\begin{align*}
    p_1'(t_{n+1}) =& \frac{6(\Delta t)^2 + 3(\Delta t)^2 + 2(\Delta t)^2}{6(\Delta t)^3} = \frac{11}{6\Delta t} \\
  p_2'(t_{n+1}) =& -\frac{0(\Delta t)^2 + 0(\Delta t)^2 + 6(\Delta t)^2}{2(\Delta t)^3} = -\frac{6}{2\Delta t} \\
  p_3'(t_{n+1}) =& \frac{3(\Delta t)^2 + 0(\Delta t)^2 + 0(\Delta t)^2}{2 (\Delta t)^3} = \frac{3}{2 \Delta t} \\
  p_4'(t_{n+1}) =& -\frac{0(\Delta t)^2 + 0(\Delta t)^2 + 2(\Delta t)^2}{6 (\Delta t)^3} = -\frac{1}{3\Delta t}
\end{align*}
From this we can construct a method to approximate $\bvec{y}_{n+1}$ which we equate with $f(t_{n+1},\bvec{y}_{n+1}):$
\begin{equation*}
P'(t_{n+1}) = \bvec{y}_{n+1}\frac{11}{6\Delta t} - \bvec{y}_{n}\frac{3}{\Delta t} + \bvec{y}_{n-1}\frac{3}{2\Delta t} - \bvec{y}_{n-2}\frac{1}{3\Delta t} = f(t_{n+1},\bvec{y}_{n+1})
\end{equation*}
We divide out $11/(6\Delta t)$ to isolate the $y_{n+1}$-term and obtain the final form of the BDF3-formula:
\begin{equation*}
  \bvec{y}_{n+1} - \bvec{y}_{n}\frac{18}{11} + \bvec{y}_{n-1}\frac{ 9}{11} - \bvec{y}_{n-2}\frac{2}{11} = \frac{6\Delta t}{11}f(t_{n+1},\bvec{y}_{n+1})
\end{equation*}
While this provides a great way to obtain higher order accurate solutions at the expense of only storage, not computational effort, we are tied to a fixed step size of $\Delta t$ for any of the formulas to be consistent.

However, sometimes we might want to change the time step size because it is much smaller than the relevant scale of the solution. One way to alleviate expanding too much computational effort is to scale up or down the order of the BDF formula. This saves some computations, but the real bottle-neck will always come from having to solve the non-linear system involving $\bvec{y}_{n+1} - c_{n+1}\Delta t f(t_{n+1}, \bvec{y}_{n+1}) + \sum_{i=0}^{N}c_{n-i} \bvec{y}_{n-i}=0.$
Fewer solutions to this system can only be made possible by making $\Delta t$ larger, so ideally we would like a way to change the time step size without ruining the complete method.
One idea is the following. When we reach a stage where we want to change the time step size and have memory of the points $\bvec{y}_{n+1},\bvec{y}_n,\bvec{y}_{n-1}$ and $\bvec{y}_{n-2},$ we construct a fourth-degree interpolation polynomial as we did before:
\begin{align*}
  P(t) = \bvec{y}_{n+1}p_1(t) + \bvec{y}_{n}p_2(t) + \bvec{y}_{n-1}p_3(t) + \bvec{y}_{n-2}p_4(t).
\end{align*}
We then determine time step size we want based on some criteria. We probably should cap it to some maximum increase, e.g. $\Delta t^{new} \leq 1.2 \Delta t$. This creates a new grid of time points:
\begin{equation*}
  t_{n+1}^{new} = t_{n+1}, \quad t_n^{new} = t_{n+1}-\Delta t^{new}, \quad t_{n-1}^{new} = t_{n+1}- 2\Delta t^{new},\quad t_{n-2}^{new} = t_{n+1}- 3\Delta t^{new}.
\end{equation*}
We can evaluate $P(t)$ on each of these points to find new approximations to $\bvec{y}_n^{new}=P(t_n^{new}), \bvec{y}_{n-1}=P(t_{n-1}^{new})$ and $\bvec{y}_{n-2}=P(t_{n-2}^{new}).$
Since the interpolation polynomial is of degree 3, our approximation will be third-order consistent, just like our BDF-3 method.

\subsection{consistency of BDF method}
If $f(t, y) = -\lambda y,$ and our first 3 solutions are exact, then we have e.g.
\begin{equation*}
  y_{n} = \exp(-\lambda 4 \Delta t), \quad   y_{n-1} = \exp(-\lambda 3 \Delta t), \quad   y_{n-2} = \exp(-\lambda 2 \Delta t).
\end{equation*}
Clearly the next point should be $\exp(-\lambda 5 \Delta t).$ The BDF method we discovered previously gives us

\begin{equation*}
  y_{n+1} - \frac{18}{11}\exp(-4\lambda \Delta t) + \frac{ 9}{11}\exp(-3\lambda \Delta t) - \frac{2}{11}\exp(-2\lambda \Delta t) = -\frac{6\Delta t}{11}\lambda y_{n+1}
\end{equation*}
Isolating for $y_{n+1}$ gives us
\begin{equation*}
  y_{n+1} = \frac{\frac{18}{11}\exp(-4\lambda \Delta t) - \frac{ 9}{11}\exp(-3\lambda \Delta t) + \frac{2}{11}\exp(-2\lambda \Delta t)}{1 + \frac{6 \lambda\Delta t}{11}}
\end{equation*}
To analyze the consistency further, we need to Taylor-expand the exponentials:
\begin{equation*}
  \exp(-n\lambda \Delta t) \approx 1 - n\lambda \Delta t + \frac{1}{2}(n\lambda \Delta t)^2 - \frac{1}{6}(n\lambda \Delta t)^3 + \frac{1}{24}(n\lambda \Delta t)^4 + \mathcal{O}(\Delta t)^5
\end{equation*}
Tallying all the contributions leads to
\begin{align*}
  y_{n+1} = (1+\frac{6\lambda \Delta t}{11})^{-1}&\left[+\frac{18}{11}\left(1 - 4\lambda \Delta t + 8(\lambda \Delta t)^2 - \frac{64}{6}(\lambda \Delta t)^3 + \frac{256}{24}(\lambda \Delta t)^4)\right)\right. \\
                                                 &\left.-\frac{9}{11}\left(1 - 3\lambda \Delta t + \frac{9}{2}(\lambda \Delta t)^2 - \frac{27}{6}(\lambda \Delta t)^3 + \frac{81}{24}(\lambda \Delta t)^4\right) \right. \\
  &+\left.\frac{2}{11}\left( 1 - 2\lambda \Delta t + 2(\lambda \Delta t)^2 - \frac{8}{6}(\lambda \Delta t)^3 + \frac{16}{24}(\lambda \Delta t)^4 \right)\right] + \mathcal{O}(\Delta t)^5
\end{align*}
We first collect all terms of equal order on the RHS inside square brackets:
\begin{align*}
  &(18 - 9 + 2)/11 = 11/11 = 1\\
  &(-72 + 27 - 4)/11 = -49/11 \\
  &(144 - 81/2 + 4) = 107.5/11 = 215/22\\
  &(-1152+243-16)/66 = -925/66 \\
  &(4608 - 729 + 32)/264 = 3911/264
\end{align*}
We also need to approximate $(1+6\lambda \Delta t/11)^{-1}$ up to fifth order:
\begin{equation*}
  \left(1+\frac{6}{11}\lambda \Delta t\right)^{-1} = 1 - \frac{6}{11}\lambda \Delta t + \left(\frac{6}{11}\lambda \Delta t\right)^2 - \left(\frac{6}{11}\lambda \Delta t\right)^3 + \left(\frac{6}{11}\lambda \Delta t\right)^4 + \mathcal{O}(\Delta t)^5
\end{equation*}
Multiplying these two leads to
\begin{align*}
  &\left[1 - \frac{6}{11}\lambda \Delta t + \left(\frac{6}{11}\lambda \Delta t\right)^2 - \left(\frac{6}{11}\lambda \Delta t\right)^3 + \left(\frac{6}{11}\lambda \Delta t\right)^4 + \mathcal{O}(\Delta t)^5\right]\cdot \\
  &\left[1 - \frac{49}{11}\lambda \Delta t + \frac{215}{22}(\lambda \Delta t)^2 - \frac{925}{66}(\lambda \Delta t)^3 + \frac{3911}{264}(\lambda \Delta t)^4 + \mathcal{O}(\Delta t)^5\right] \\
  &= 1 + (\lambda \Delta t)\left(\frac{-6-49}{11}\right) + (\lambda \Delta t)^2\left(\frac{6^2}{11^2} + \frac{215}{22} + \frac{294}{11^2}\right) \\
  &+ (\lambda \Delta t)^3\left(-\frac{216}{11^3} - \frac{36\cdot 49}{11^3} - \frac{215\cdot 6}{22\cdot 11} - \frac{925}{66}\right) \\
  &+ (\lambda \Delta t)^4\left(\frac{6^4}{11^4} + \frac{49\cdot 6^3}{11^4} + \frac{6^2\cdot 215}{11^2\cdot 22}+\frac{925\cdot 6}{66\cdot 11} + \frac{3911}{264}\right) + \mathcal{O}(\Delta t)^5
\end{align*}
Simplifying the terms leads to
\begin{align*}
  &= 1 - (\lambda \Delta t)\left(\frac{55}{11}\right) + (\lambda \Delta t)^2\left(\frac{3025}{242} \right) - (\lambda \Delta t)^3\left(\frac{166375}{66\cdot 11^2} \right) \\
  &+ (\lambda \Delta t)^4\left( \frac{55191246}{2108304}\right) + \mathcal{O}(\Delta t)^5
\end{align*}
We can simplify the fractions down to
\begin{align*}
  &= 1 - 5(\lambda \Delta t) + \frac{1}{2}(5\lambda \Delta t)^2 - \frac{1}{6}(5 \lambda \Delta t)^3  + \mathcal{O}(\Delta t)^4
\end{align*}
The first four terms match the Taylor series of $\exp(-5\lambda \Delta t)$ but the fifth term is off. We expect $5^4/24$ for the final term but have instead $6911/(11\cdot 24).$ In other words, we obtained the approximation $5^4 = 625 \approx 628.27 = 6911/11.$
This means our BDF-3 formula is third-order consistent.

\subsection{Solving BDF equations}
To solve a system like the BDF equations, we perform the following steps. Consider again the BDF-3 method:
\begin{equation*}
  \bvec{y}_{n+1} - \frac{18}{11}\bvec{y}_n + \frac{9}{11}\bvec{y}_{n-1} - \frac{2}{11}\bvec{y}_{n-2} = \frac{6 \Delta t}{11}f(t_{n+1}, \bvec{y}_{n+1})
\end{equation*}
which we can rewrite as
\begin{equation*}
  \bvec{y}_{n+1} - \frac{18}{11}\bvec{y}_n + \frac{9}{11}\bvec{y}_{n-1} - \frac{2}{11}\bvec{y}_{n-2} - \frac{6 \Delta t}{11}\bvec{f}(t_{n+1}, \bvec{y}_{n+1}) = \bvec{0}.
\end{equation*}
We can straightforwardly apply Newton's method to this since we see that
\begin{equation*}
  \bvec{I}  - \frac{6 \Delta t}{11}\bvec{J}(t_{n+1}, \bvec{y}_{n+1})
\end{equation*}
is the full Jacobi matrix of the BDF method, where $\bvec{J}$ is the Jacobi matrix of the ODE. The Newton type iterations therefore become
\begin{equation*}
  \bvec{y}_{n+1}^{k+1} = \bvec{y}_{n+1}^k + \left(\bvec{I}  - \frac{6 \Delta t}{11}\bvec{J}(t_{n+1}, \bvec{y}^{k}_{n+1})\right)^{-1}\left( \bvec{y}^{k}_{n+1} - \frac{18}{11}\bvec{y}_n + \frac{9}{11}\bvec{y}_{n-1} - \frac{2}{11}\bvec{y}_{n-2} - \frac{6 \Delta t}{11}\bvec{f}(t_{n+1}, \bvec{y}_{n+1}^{k}) \right)
\end{equation*}












\section{General linear methods}
A generalization to Runge-Kutta methods is the family of so-called general linear methods (GLM).
They can be thought of as a combination of RK and multistep methods.
A GLM consists of $s$ stages and $r$ history points. For $r=1$ we recover normal Runge-Kutta methods.
The $s$ stage points $Y_i,~0<i<s-1$ and stage derivatives $F_i$ are defined as follows:
\begin{align*}
  Y_i = \sum_{j=0}^{s-1} a_{ij}\Delta t F_j + \sum_{j=0}^{r-1}u_{ij} y_{n-j}, \qquad F_j := f(t+\Delta t c_i, Y_j).
\end{align*}
As one can see, we now have two matrices $\bvec{A} = (a_{ij})$ and $\bvec{U} = (u_{ij})$ that define the stages.
For the update to the numerical solution we similarly have two vectors $\bvec{b}$ and $\bvec{v}$:
\begin{align*}
  y_{n+1} = \sum_{i=0}^{s-1} b_i\Delta t F_i + \sum_{i=0}^{r-1}v_i y_{n-i}.
\end{align*}
Applying a GLM to the test problem $f(t,y) = \lambda y$ leads to the following:
\begin{align*}
  Y_i =& \sum_{j=0}^{s-1} z a_{ij} Y_j + \sum_{j=0}^{r-1} u_{ij}y_{n-j} \\
  y_{n+1} =& \sum_{i=0}^{s-1} z b_i Y_i + \sum_{i=0}^{r-1} v_iy_{n-i},
\end{align*}
where $z:=\lambda \Delta t.$
Applying some linear algebra we can write the stage values the following way:
\begin{align*}
  \bvec{Y} := \begin{pmatrix}
    Y_0 \\
    Y_1 \\
    Y_2 \\
    \vdots \\
    Y_{s-1}
  \end{pmatrix} = z \bvec{A} \bvec{Y} + \bvec{U}\begin{pmatrix}
    y_n \\
    y_{n-1} \\
    y_{n-2} \\
    \vdots \\
    y_{n-r+1}
  \end{pmatrix},
\end{align*}
which can be compacted into $(\bvec{I} - z \bvec{A})\bvec{Y} = \bvec{U} \bvec{y},$ with $\bvec{y} = (y_n, y_{n-1}, y_{n-2}, \hdots, y_{n-r+1})^T.$
Thus, for the stage points we simply find
\begin{equation*}
  \bvec{Y} = (\bvec{I} - z\bvec{A})^{-1}\bvec{Uy}
\end{equation*}
and hence the update becomes
\begin{equation*}
  \bvec{y}(t_{n+1}) = z \bvec{b} \cdot (\bvec{I} - z\bvec{A})^{-1}\bvec{Uy} + \bvec{v}\cdot\bvec{y}.
\end{equation*}
As a check, consider $\bvec{U}=1$ and $\bvec{v} = 1$ and $r=1,$ in which case we should recover the result for RK methods.
We see that in this case
\begin{equation*}
  \bvec{y}(t_{n+1}) = z \bvec{b} \cdot (\bvec{I} - z\bvec{A})^{-1}\bvec{y}(t_n) + \bvec{y}(t_n).
\end{equation*}
If we have $y(t_n) = 1$ and divide by it, we recover the stability function:
\begin{equation*}
  y(t_{n+1})/y(t_n) = 1 + z \bvec{b} \cdot (\bvec{I} - z\bvec{A})^{-1}(1,1,1,\hdots,1)^T,
\end{equation*}
which indeed is the stability function for RK methods.






%\section{Optimizations in implementing IRK methods}
%Consider again a three-stage IRK method, e.g. a three-stage Lobatto IIIC method.
%We have
%\begin{align*}
%  \bvec{Y} :=& \Delta t
%  \begin{pmatrix}
%    a_{11}\bvec{k}_1 + a_{12}\bvec{k}_2 + a_{13}\bvec{k}_3 \\
%    a_{21}\bvec{k}_1 + a_{22}\bvec{k}_2 + a_{23}\bvec{k}_3 \\
%    a_{31}\bvec{k}_1 + a_{32}\bvec{k}_2 + a_{33}\bvec{k}_3
%  \end{pmatrix} \\
%  =& \Delta t \begin{pmatrix}
%    a_{11} & 0 & 0 &  a_{12} & 0 & 0 & a_{13} & 0 & 0 \\
%    0 & a_{11} & 0 & 0 &  a_{12} & 0 & 0 & a_{13} & 0 \\
%    0 & 0 & a_{11} & 0 & 0 &  a_{12} & 0 & 0 & a_{13} \\
%    a_{21} & 0 & 0 &  a_{22} & 0 & 0 & a_{23} & 0 & 0 \\
%    0 & a_{21} & 0 & 0 &  a_{22} & 0 & 0 & a_{23} & 0 \\
%    0 & 0 & a_{21} & 0 & 0 &  a_{22} & 0 & 0 & a_{23} \\
%    a_{31} & 0 & 0 &  a_{32} & 0 & 0 & a_{33} & 0 & 0 \\
%    0 & a_{31} & 0 & 0 &  a_{32} & 0 & 0 & a_{33} & 0 \\
%    0 & 0 & a_{31} & 0 & 0 &  a_{32} & 0 & 0 & a_{33}
%     \end{pmatrix}
%     \begin{pmatrix}
%       k_{11} \\
%       k_{12} \\
%       k_{13} \\
%       k_{21} \\
%       k_{22} \\
%       k_{23} \\
%       k_{31} \\
%       k_{32} \\
%       k_{33} \\
%     \end{pmatrix}
%\end{align*}
%With this we solved the equation
%\begin{align*}
%  (\Delta t \bvec{A}\otimes \bvec{I})^{-1}\bvec{Y} = \begin{pmatrix}
%  \bvec{f}(t + c_1\Delta t, \bvec{y} + \bvec{Y}_1) \\
%  \bvec{f}(t + c_2\Delta t, \bvec{y} + \bvec{Y}_2) \\
%  \bvec{f}(t + c_3\Delta t, \bvec{y} + \bvec{Y}_3) \\
%  \end{pmatrix}
%\end{align*}
%which written out explicitly is
%\begin{align*}
%  \Delta t \begin{pmatrix}
%    a_{11} & 0 & 0 &  a_{12} & 0 & 0 & a_{13} & 0 & 0 \\
%    0 & a_{11} & 0 & 0 &  a_{12} & 0 & 0 & a_{13} & 0 \\
%    0 & 0 & a_{11} & 0 & 0 &  a_{12} & 0 & 0 & a_{13} \\
%    a_{21} & 0 & 0 &  a_{22} & 0 & 0 & a_{23} & 0 & 0 \\
%    0 & a_{21} & 0 & 0 &  a_{22} & 0 & 0 & a_{23} & 0 \\
%    0 & 0 & a_{21} & 0 & 0 &  a_{22} & 0 & 0 & a_{23} \\
%    a_{31} & 0 & 0 &  a_{32} & 0 & 0 & a_{33} & 0 & 0 \\
%    0 & a_{31} & 0 & 0 &  a_{32} & 0 & 0 & a_{33} & 0 \\
%    0 & 0 & a_{31} & 0 & 0 &  a_{32} & 0 & 0 & a_{33}
%     \end{pmatrix}^{-1}
%     \begin{pmatrix}
%       Y_{11} \\
%       Y_{12} \\
%       Y_{13} \\
%       Y_{21} \\
%       Y_{22} \\
%       Y_{23} \\
%       Y_{31} \\
%       Y_{32} \\
%       Y_{33} \\
%     \end{pmatrix} = \\
%  \begin{pmatrix}
%      f_1(t+c_1\Delta t, y_1+Y_{11}+Y_{12}+Y_{13}) \\
%      f_2(t+c_1\Delta t, y_2+Y_{11}+Y_{12}+Y_{13}) \\
%      f_3(t+c_1\Delta t, y_3+Y_{11}+Y_{12}+Y_{13}) \\
%      f_1(t+c_2\Delta t, y_1+Y_{21}+Y_{22}+Y_{23}) \\
%      f_2(t+c_2\Delta t, y_2+Y_{21}+Y_{22}+Y_{23}) \\
%      f_3(t+c_2\Delta t, y_3+Y_{21}+Y_{22}+Y_{23}) \\
%      f_1(t+c_3\Delta t, y_1+Y_{31}+Y_{32}+Y_{33}) \\
%      f_2(t+c_3\Delta t, y_2+Y_{31}+Y_{32}+Y_{33}) \\
%      f_3(t+c_3\Delta t, y_3+Y_{31}+Y_{32}+Y_{33}) \\
%  \end{pmatrix}
%\end{align*}
%Now, suppose we can find a transformation for the matrix $\bvec{A}$ such that
%\begin{equation*}
%  \bvec{A} = \bvec{T}^{-1}\bvec{U}\bvec{T}
%\end{equation*}


\section{Some more test equations}
Other good tests for our solver would be the systems of equations obtained after discretizing a PDE. This wll also provide a good way to make sure sparse matrix support is working.
Consider the homogeneous heat equation
\begin{equation*}
  \frac{\partial u}{\partial t} = \kappa \left[ \frac{\partial^2u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right]
\end{equation*}
subject to boundary conditions $u(x, 0) = u_0,$ $u(x, L_y)=u_0$, $u(0, y) = u_0$ and $\partial u/\partial x(L_x, y) = -\Gamma.$

We can discretize by means of finite (central) differences, but of course we can also directly compute the solution if we are given an initial value.
To find an exact solution we apply separation of variables and obtain
\begin{equation*}
  T'/T = \kappa (X''/X + Y''/Y) = const. := -\lambda
\end{equation*}
and we find $T(t) = A \exp(-\lambda t).$
Similarly, we find
\begin{equation*}
  \kappa X''/X = -\lambda - \kappa Y''/Y = const.
\end{equation*}
and so $X(x) =  A \cos(kx) + B \sin(kx)$ and $Y(y) = A \cos(ky) + B \sin(ky).$
Substituting all this back into the PDE gives us
\begin{align*}
  -\lambda &\exp(-\lambda t) \left[A\cos(kx) + B\sin(kx) + C\cos(ky) + D\sin(ky)\right] \\
  =&\kappa \exp(-\lambda t) \left[ -k_x^2 A\cos(k_xx)  -k_x^2B\cos(k_xx)\right]\left[by +  C\cos(k_yy) + D\sin(k_yy)\right] \\
           &+ \kappa \exp(-\lambda t) \left[ ax +  A\cos(k_xx)  + B\cos(k_xx)\right]\left[-k_y^2C\cos(k_yy) - k_y^2 D\sin(k_yy)\right] \\
  =& -\kappa (k_x^2  + k_y^2) u = -\lambda u
\end{align*}
For this to be equal we clearly must have $k_x^2  + k_y^2 = \lambda.$
Both $k_x$ and $k_y$ can only take on specific values which follow from the boundary conditions.
For $y=0$ we have $u(x,0) = u_0$ and for $y = L_y$ we have $u(x,L_y) = u_0.$ This means that $k_y= n \pi$ with $n \in \mathbb{R}$ and we can write
\begin{equation*}
  Y(y) = \sum_{n=0}^\infty (A_n\cos(n \pi y) + B_n \sin(n\pi y))
\end{equation*}
and $A_0 = u_0.$

\clearpage{}
\begin{appendix}
  \section{More coefficients for different methods}
  %recommended methods
  %LOBATTO_IIIC_85 (FIX)
  %rest
  %LOBATTO_IIIA_32
  %RADAU_IIA_53
  %RADAU_IIA_95
  %RADAU_IIA_137
  %GAUSS_LEGENDRE_147
  %LOBATTO_IIIA_127
  %LOBATTO_IIIC_127

  \subsection{Lobatto IIIC methods}
  For the Lobatto IIIC methods of 4 and 5 stages, we can derive similar embedded pairs.
  The 4 stage method has Butcher tableau
  \begin{table}[h]
    \centering
    \caption{4-stage Lobatto IIIC method}
      \def\arraystretch{1.5}%
  \begin{tabular}{c|cccc}
    0                             & $\frac{1}{12}$ & $-\frac{\sqrt{5}}{12}$ & $\frac{\sqrt{5}}{12}$ & $-\frac{1}{12}$ \\
    $\half - \frac{\sqrt{5}}{10}$ & $\frac{1}{12}$ & $\frac{1}{4}$ & $\frac{10-7\sqrt{5}}{60}$ & $\frac{\sqrt{5}}{60} $ \\
    $\half + \frac{\sqrt{5}}{10}$ & $\frac{1}{12}$ & $\frac{10+7\sqrt{5}}{60}$ & $\frac{1}{4}$ & $-\frac{\sqrt{5}}{60}$ \\
    1                             & $\frac{1}{12}$ & $\frac{5}{12}$ & $\frac{5}{12}$ & $\frac{1}{12}$ \\ \hline
    {} & $\frac{1}{12}$ & $\frac{5}{12}$ & $\frac{5}{12}$ & $\frac{1}{12}$ \\
    {} & $\frac{1}{12}$ & $\frac{1}{12}$ & $\frac{5}{12}$ & $\frac{5}{12}$ \\
  \end{tabular}
\end{table}

For the 5-stage, we apply the simple cyclic rotation (each $\bvec{b}$ coefficient shifts one to the right). This leads to a method that is $A(\alpha)$ stable with $\alpha\approx 71\degree,$ which is not spectacular
  \begin{table}[h]
    \centering
    \caption{5-stage Lobatto IIIC method}
    \def\arraystretch{1.5}%
  \begin{tabular}{c|ccccc}
    0                             & $\frac{1}{20}$ & $-\frac{7}{60}$ & $\frac{2}{15}$ & $-\frac{7}{60}$ & $\frac{1}{20}$ \\
    $\half - \frac{\sqrt{21}}{14}$ & $\frac{1}{20}$ & $\frac{29}{180}$ & $\frac{47-15\sqrt{21}}{315}$ & $\frac{203-30\sqrt{21}}{1260}$ & $-\frac{3}{140}$ \\
    $\half$ & $\frac{1}{20}$ & $\frac{329+105\sqrt{21}}{2880}$ & $\frac{73}{360}$ & $\frac{329-105\sqrt{21}}{2880}$ & $\frac{3}{160}$ \\
    $\half + \frac{\sqrt{21}}{14}$ & $\frac{1}{20}$ & $\frac{203+30\sqrt{21}}{1260}$ & $\frac{47+15\sqrt{21}}{315}$ & $\frac{29}{180} $ & $-\frac{3}{140}$ \\
    $1$ & $\frac{1}{20}$ & $\frac{49}{180}$ & $\frac{16}{45}$ & $\frac{49}{180}$ & $\frac{1}{20}$ \\ \hline
    {} & $\frac{1}{20}$ & $\frac{49}{180}$ & $\frac{16}{45}$ & $\frac{49}{180}$ & $\frac{1}{20}$ \\
    {} & $\frac{1}{20}$ & $\frac{1}{20}$ & $\frac{49}{180}$ & $\frac{16}{45}$ & $\frac{49}{180}$
  \end{tabular}
\end{table}

  \subsection{Radau IIA methods}
\end{appendix}



\end{document}