From 500472e68826f065006ea8645970413d93aca7b5 Mon Sep 17 00:00:00 2001
From: Andrey A Popov <apopov@vt.edu>
Date: Tue, 29 Aug 2023 21:45:17 -0500
Subject: [PATCH 01/20] First draft of lorenz 96 docs

---
 docs/references.bib                          | 15 +++++-
 toolbox/+otp/+lorenz96/+presets/Canonical.m  | 23 ++++++---
 toolbox/+otp/+lorenz96/+presets/PopovSandu.m | 24 ++++++++-
 toolbox/+otp/+lorenz96/Lorenz96Parameters.m  |  6 +--
 toolbox/+otp/+lorenz96/Lorenz96Problem.m     | 51 ++++++++++++++++++--
 5 files changed, 102 insertions(+), 17 deletions(-)

diff --git a/docs/references.bib b/docs/references.bib
index 141c7a39..af3d0a2f 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -105,5 +105,18 @@ @inproceedings{Lor96
   volume={1},
   number={1},
   year={1996},
-  organization={Reading}
+  organization={Reading},
+  doi={10.1017/CBO9780511617652.004}
+}
+
+@article{PS19,
+  title={A Bayesian approach to multivariate adaptive localization in ensemble-based data assimilation with time-dependent extensions},
+  author={Popov, Andrey A and Sandu, Adrian},
+  journal={Nonlinear Processes in Geophysics},
+  volume={26},
+  number={2},
+  pages={109--122},
+  year={2019},
+  publisher={Copernicus Publications G{\"o}ttingen, Germany},
+  doi={10.5194/npg-26-109-2019}
 }
diff --git a/toolbox/+otp/+lorenz96/+presets/Canonical.m b/toolbox/+otp/+lorenz96/+presets/Canonical.m
index 38c0b0fc..ac4c3f56 100644
--- a/toolbox/+otp/+lorenz96/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz96/+presets/Canonical.m
@@ -1,12 +1,23 @@
 classdef Canonical < otp.lorenz96.Lorenz96Problem
-    %CANONICAL This is the original problem presented in the literature.
-    %  The initial condition is a slight perturbation of a critical point.
-    %
-    % See:
-    %     Lorenz, E. N. (1996, September). Predictability: A problem partly solved. 
-    %     In Proc. Seminar on predictability (Vol. 1, No. 1).
+    % Original Lorenz '96 preset presented in :cite:p:`Lor96`
+    % which uses time span $t \in [0, 720]$, $N = 40$, $F=8$, and initial
+    % conditions of $y_i = 8$ for all $i$ except for $y_20=8.008$. 
+
     methods
         function obj = Canonical(varargin)
+            % Create a Canonial problem object.
+            %
+            % Parameters
+            % ----------
+            % Size : numeric(1, 1)
+            %    The size of the problem as a positive integer.
+            % Forcing : numeric(N, 1)
+            %    The forcing as a vector of N constants.
+            %
+            % Returns
+            % -------
+            % obj : Canonial
+            %    The constructed problem.
 
             p = inputParser;
             p.addParameter('Size', 40, @isscalar);
diff --git a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
index f04fdada..3c171f1d 100644
--- a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
+++ b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
@@ -1,9 +1,29 @@
 classdef PopovSandu < otp.lorenz96.Lorenz96Problem
-    %POPOVSANDU
-    % Used in https://doi.org/10.5194/npg-26-109-2019
+    % A preset that has a cyclic forcing function that is different for
+    % every variable. This preset was created for :cite:p:`PS19`.
+    % This preset uses time span $t \in [0, 720]$, $N = 40$, $F=8$, 
+    % four partitions, a forcing period of one time unit, and initial
+    % conditions of $y_i = 8$ for all $i$ except for $y_20=8.008$. 
     
     methods
         function obj = PopovSandu(varargin)
+            % Create a PopovSandu problem object.
+            %
+            % Parameters
+            % ----------
+            % Size : numeric(1, 1)
+            %    The size of the problem as a positive integer.
+            % Partitions : numeric(1, 1)
+            %    The number of partitions into which to divide the
+            %    variables.
+            % ForcingPeriod : numeric(1, 1)
+            %    The period of the forcing function in radians per unit
+            %    time.
+            %
+            % Returns
+            % -------
+            % obj : PopovSandu
+            %    The constructed problem.
             
             p = inputParser;
             
diff --git a/toolbox/+otp/+lorenz96/Lorenz96Parameters.m b/toolbox/+otp/+lorenz96/Lorenz96Parameters.m
index c1e94338..bb3ee4a6 100644
--- a/toolbox/+otp/+lorenz96/Lorenz96Parameters.m
+++ b/toolbox/+otp/+lorenz96/Lorenz96Parameters.m
@@ -1,10 +1,8 @@
 classdef Lorenz96Parameters
-    %LORENZ96PARAMETERS User-configurable parameters for the Lorenz 96 problem
-    %
-    %   See also otp.lorenz96.Lorenz96Problem
+    % Parameters for the Lorenz '96 problem.
     
     properties
-        %F is a forcing function, scalar, or vector
+        % A forcing function, scalar, or vector.
         F %MATLAB ONLY: {mustBeA(F, {'numeric','function_handle'})}
     end
 end
diff --git a/toolbox/+otp/+lorenz96/Lorenz96Problem.m b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
index b8e26883..72b3fd93 100644
--- a/toolbox/+otp/+lorenz96/Lorenz96Problem.m
+++ b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
@@ -1,7 +1,35 @@
 classdef Lorenz96Problem <  otp.Problem
-    %LORENZ96PROBLEM  The Lorenz 96 model is a classic model for testing data assimilation techniques.
-    % 
-    % See also otp.loren96.Lorenz96Parameters
+    % A chaotic system modeling the transfer of some quantity along some
+    % longitude.
+    %
+    %
+    % The $N$ variable dynamics :cite:p:`Lor96` are represented by the equation,
+    %
+    % $$
+    % y_i' = -y_{i-1}\left(y_{i-2} - y_{i+1}\right) - y_i + f(t),\quad i \in \mathbb{Z}_N,
+    % $$
+    %
+    % that exhibits chaotic behavior for certain values of the forcing function $f$.
+    %
+    % Notes
+    % -----
+    % +---------------------+-----------------------------------------------------------+
+    % | Type                | ODE                                                       |
+    % +---------------------+-----------------------------------------------------------+
+    % | Number of Variables | $N$                                                       |
+    % +---------------------+-----------------------------------------------------------+
+    % | Stiff               | no                                                        |
+    % +---------------------+-----------------------------------------------------------+
+    %
+    % Example
+    % -------
+    % >>> problem = otp.lorenz96.presets.Canonical('Forcing', @(t) 8 + 4*sin(t));
+    % >>> sol = problem.solve();
+    % >>> problem.plotPhaseSpace(sol);
+    %
+    % See also
+    % --------
+    % :doc:`Lorenz '63 Problem </problems/Lorenz63Problem>`
     
     properties (SetAccess = private)
         DistanceFunction
@@ -9,6 +37,21 @@
     
     methods
         function obj = Lorenz96Problem(timeSpan, y0, parameters)
+            % Create a Lorenz '96 problem object.
+            %
+            % Parameters
+            % ----------
+            % timeSpan : numeric(1, 2)
+            %    The start and final time.
+            % y0 : numeric(N, 1)
+            %    The initial conditions.
+            % parameters : Lorenz96Parameters
+            %    The parameters.
+            %
+            % Returns
+            % -------
+            % obj : Lorenz96Problem
+            %    The constructed problem.
             obj@otp.Problem('Lorenz 96', [], timeSpan, y0, parameters);
         end
     end
@@ -37,7 +80,7 @@
                 'Vectorized', 'on');
             
             % We also provide a canonical distance function as is standard for
-            % localisation in Data Assimilation. This is heavily tied to this
+            % localization in Data Assimilation. This is heavily tied to this
             % problem.
             
             obj.DistanceFunction = @(t, y, i, j) otp.lorenz96.distanceFunction(t, y, i, j);

From cd53d2bc3f82a7cc6a97f8206d00bcb7d498612c Mon Sep 17 00:00:00 2001
From: Arash Sarshar <elswit@users.noreply.github.com>
Date: Tue, 29 Aug 2023 19:55:01 -0700
Subject: [PATCH 02/20] started ascher DAE doc writeup

---
 docs/references.bib                           | 14 ++++
 notebooks/quick-start.ipynb                   | 62 ++++++++++++++++--
 .../+ascherlineardae/AscherLinearDAEProblem.m | 64 ++++++++++++++++++-
 3 files changed, 134 insertions(+), 6 deletions(-)

diff --git a/docs/references.bib b/docs/references.bib
index 141c7a39..c2312d9d 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -107,3 +107,17 @@ @inproceedings{Lor96
   year={1996},
   organization={Reading}
 }
+
+
+@article{Asc89,
+author = {Ascher, Uri},
+title = {On Symmetric Schemes and Differential-Algebraic Equations},
+journal = {SIAM Journal on Scientific and Statistical Computing},
+volume = {10},
+number = {5},
+pages = {937-949},
+year = {1989},
+doi = {10.1137/0910054}
+}
+
+
diff --git a/notebooks/quick-start.ipynb b/notebooks/quick-start.ipynb
index 00b490aa..73ac5a2b 100644
--- a/notebooks/quick-start.ipynb
+++ b/notebooks/quick-start.ipynb
@@ -181,6 +181,37 @@
     "The parameters of a model are stored as properties in the problem object:"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\u001b[?2004l\n",
+      "problem =\n",
+      "\n",
+      "  otp.lorenz63.presets.Canonical object with properties:\n",
+      "\n",
+      "            Name: Lorenz Equations\n",
+      "         NumVars: [1x1 double]\n",
+      "      Parameters: [1x1 otp.lorenz63.Lorenz63Parameters]\n",
+      "             RHS: [1x1 otp.RHS]\n",
+      "        TimeSpan: [1x2 double]\n",
+      "              Y0: [3x1 double]\n",
+      "\n",
+      "\u001b[?2004h"
+     ]
+    }
+   ],
+   "source": [
+    "problem"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 4,
@@ -528,7 +559,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 45,
+   "execution_count": 5,
    "metadata": {
     "tags": []
    },
@@ -552,7 +583,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 51,
+   "execution_count": 8,
    "metadata": {
     "tags": []
    },
@@ -569,6 +600,7 @@
       "\u001b[?2004h\u001b[?2004l\n",
       "\u001b[?2004h\u001b[?2004l\n",
       "\u001b[?2004h\u001b[?2004l\n",
+      "\u001b[?2004h\u001b[?2004l\n",
       "\u001b[?2004h"
      ]
     },
@@ -591,9 +623,31 @@
     "\n",
     "\n",
     "subplot(132); contourf(reshape(sol.y(:,end), dim_x, []), 10)\n",
-    "axis('equal'); title(sprintf( 't = %1.1f' , sol.x(end)));  xlabel('x'); ylabel('y')\n"
+    "axis('equal'); title(sprintf( 't = %1.1f' , sol.x(end)));  xlabel('x'); ylabel('y')\n",
+    "\n",
+    "\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\u001b[?2004l\r",
+      "\r\n",
+      "error: print: no figure to print\r\n",
+      "error: called from\r\n",
+      "    print at line 463 column 5\r\n",
+      "\u001b[?2004h"
+     ]
+    }
+   ],
+   "source": []
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -603,7 +657,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 56,
+   "execution_count": 13,
    "metadata": {
     "tags": []
    },
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
index 34764a34..5d95208b 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
@@ -1,8 +1,69 @@
 classdef AscherLinearDAEProblem < otp.Problem
-    %ASCHERLINEARDAEPROBLEM This is an Index-1 DAE problem
+    % A simple linear differential algebraic problem.
+    %
+    % The Ascher linear DAE Problem :cite:p:`Asc89` is given by $M(t) y' = A(t) y + q(t) $ with
+    %
+    % $$
+    % M=\left(\begin{array}{cc}
+    % 1 & -t \\
+    % 0 & 0
+    % \end{array}\right), \quad A=\left(\begin{array}{cc}
+    % -1 & 1+t \\
+    % \beta & -1-\beta t
+    % \end{array}\right), \quad {q}=\left(\begin{array}{c}
+    % 0 \\
+    % \sin t
+    % \end{array}\right), 
+    % $$
+    %
+    % with $t \in [0,1]$, and $y_0 = [1, \beta]^T$.  The exact solution
+    % is given by 
+    %
+    % $$
+    % y = \begin{pmatrix}
+    % t \sin(t) + (1 + \beta  t) e^{-t}\\
+    % \beta  e^{-t} + \sin(t)
+    % \end{pmatrix}.
+    % $$
+    %
+    % Due to its stiffness and time-dependant mass
+    % matrix and righ-hand-side function, this simple DAE problem can 
+    % become challenging to solve. This problem is used in :cite:p:`Asc89` 
+    % to study convergence of implcit solvers. 
+    %
+    % Notes
+    % -----
+    % +---------------------+-----------------------------------------+
+    % | Type                | DAE                                     |
+    % +---------------------+-----------------------------------------+
+    % | Number of Variables | 2                                       |
+    % +---------------------+-----------------------------------------+
+    % | Stiff               | typically, depending on $\beta$         |
+    % +---------------------+-----------------------------------------+
+    %
+    % Example
+    % -------
+    % >>> problem = otp.ascherlineardae.presets.Canonical;
+    % >>> sol = problem.solve();
+    % >>> problem.plotPhaseSpace(sol);
     %
     methods
         function obj = AscherLinearDAEProblem(timeSpan, y0, parameters)
+            % Create an Ascher linear DAE problem object.
+            %
+            % Parameters
+            % ----------
+            % timeSpan : numeric(1, 2)
+            %    The start and final time.
+            % y0 : numeric(2, 1)
+            %    The initial conditions.
+            % parameters : BrusselatorParameters
+            %    The parameters.
+            %
+            % Returns
+            % -------
+            % obj : BrusselatorProblem
+            %    The constructed problem.
             obj@otp.Problem('Ascher Linear DAE', 2, timeSpan, y0, parameters);
         end
     end
@@ -35,4 +96,3 @@ function onSettingsChanged(obj)
         end
     end
 end
-

From 3d932af1229bbf2536ef639bc49c2e6f47ed99ff Mon Sep 17 00:00:00 2001
From: Arash Sarshar <elswit@users.noreply.github.com>
Date: Tue, 5 Sep 2023 17:28:08 -0700
Subject: [PATCH 03/20] docs for Ascher linear DAE

---
 .../+ascherlineardae/+presets/Canonical.m     |  9 +++----
 .../AscherLinearDAEParameters.m               |  5 ++--
 .../+ascherlineardae/AscherLinearDAEProblem.m | 26 ++++++++++++-------
 3 files changed, 22 insertions(+), 18 deletions(-)

diff --git a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
index e7ec574d..3f3432d0 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
@@ -1,10 +1,7 @@
 classdef Canonical < otp.ascherlineardae.AscherLinearDAEProblem
-    %CANONICAL The problem formulation of the linear DAE from the literature
-    %
-    % See
-    %    Ascher, Uri. "On symmetric schemes and differential-algebraic equations."
-    %    SIAM journal on scientific and statistical computing 10.5 (1989): 937-949.
-    
+    % The original problem defined by Uri Ascher in :cite:p:`Asc89` with 
+    % $\beta = 0.5$
+    % 
     methods
         function obj = Canonical(beta)
             tspan = [0.0; 1];
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
index ae660c65..b50416c0 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
@@ -1,7 +1,8 @@
 classdef AscherLinearDAEParameters
-    %ASCHERLINEARDAEPARAMETERS 
+    % Parameters for Ascher Linear DAE problem 
     properties
-        %Beta is an arbitrary scalar parameter
+        % Beta is a scalar parameter in the linear model. It also affects the initial condition
+        % of the algebraic variable.
         Beta %MATLAB ONLY: (1,1) {mustBeNumeric} = 0.5
     end
 end
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
index 5d95208b..195a709c 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
@@ -1,7 +1,8 @@
 classdef AscherLinearDAEProblem < otp.Problem
     % A simple linear differential algebraic problem.
     %
-    % The Ascher linear DAE Problem :cite:p:`Asc89` is given by $M(t) y' = A(t) y + q(t) $ with
+    % The Ascher linear DAE Problem :cite:p:`Asc89` is an index-1 differential agebraic 
+    % equation given by $M(t) y' = A(t) y + q(t) $ where
     %
     % $$
     % M=\left(\begin{array}{cc}
@@ -16,7 +17,7 @@
     % \end{array}\right), 
     % $$
     %
-    % with $t \in [0,1]$, and $y_0 = [1, \beta]^T$.  The exact solution
+    % defined on timespan $t \in [0,1]$, and initial condition $y_0 = [1, \beta]^T$. The exact solution
     % is given by 
     %
     % $$
@@ -27,9 +28,9 @@
     % $$
     %
     % Due to its stiffness and time-dependant mass
-    % matrix and righ-hand-side function, this simple DAE problem can 
-    % become challenging to solve. This problem is used in :cite:p:`Asc89` 
-    % to study convergence of implcit solvers. 
+    % matrix, this simple DAE problem can 
+    % become challenging to solve. This problem is introduced in :cite:p:`Asc89` 
+    % to study the convergence of implcit solvers applied to DAEs.
     %
     % Notes
     % -----
@@ -43,10 +44,15 @@
     %
     % Example
     % -------
-    % >>> problem = otp.ascherlineardae.presets.Canonical;
-    % >>> sol = problem.solve();
+    % >>> problem = otp.ascherlineardae.presets.Canonical(0.1);
+    % >>> sol = problem.solve('MaxStep',1e-5);
     % >>> problem.plotPhaseSpace(sol);
-    %
+
+    properties (Access = private, Constant)
+        NumComps = 2
+        VarNames = 'yz'
+    end
+
     methods
         function obj = AscherLinearDAEProblem(timeSpan, y0, parameters)
             % Create an Ascher linear DAE problem object.
@@ -57,12 +63,12 @@
             %    The start and final time.
             % y0 : numeric(2, 1)
             %    The initial conditions.
-            % parameters : BrusselatorParameters
+            % parameters : AscherLinearDAEParameters
             %    The parameters.
             %
             % Returns
             % -------
-            % obj : BrusselatorProblem
+            % obj : AscherLinearDAEProblem
             %    The constructed problem.
             obj@otp.Problem('Ascher Linear DAE', 2, timeSpan, y0, parameters);
         end

From 427216dadc813a2240452960303fa066c9e9cfba Mon Sep 17 00:00:00 2001
From: Arash Sarshar <elswit@users.noreply.github.com>
Date: Sat, 16 Sep 2023 16:02:39 -0700
Subject: [PATCH 04/20] up Ascher docs, added presets

---
 docs/references.bib                           | 11 +++++
 .../+ascherlineardae/+presets/Canonical.m     | 29 +++++++-----
 .../+otp/+ascherlineardae/+presets/Petzold.m  | 19 ++++++++
 .../+otp/+ascherlineardae/+presets/Stiff.m    | 19 ++++++++
 .../AscherLinearDAEParameters.m               |  6 +--
 .../+ascherlineardae/AscherLinearDAEProblem.m | 45 +++++++++----------
 6 files changed, 92 insertions(+), 37 deletions(-)
 create mode 100644 toolbox/+otp/+ascherlineardae/+presets/Petzold.m
 create mode 100644 toolbox/+otp/+ascherlineardae/+presets/Stiff.m

diff --git a/docs/references.bib b/docs/references.bib
index c2312d9d..be5d7598 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -120,4 +120,15 @@ @article{Asc89
 doi = {10.1137/0910054}
 }
 
+@article{Pet86,
+ author = {L. R. Petzold},
+ journal = {SIAM Journal on Numerical Analysis},
+ number = {4},
+ pages = {837--852},
+ publisher = {Society for Industrial and Applied Mathematics},
+ title = {Order Results for Implicit Runge-Kutta Methods Applied to Differential/ Algebraic Systems},
+ volume = {23},
+ year = {1986},
+ url = {https://www.jstor.org/stable/2157625}
+}
 
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
index 3f3432d0..865a3c94 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
@@ -1,19 +1,28 @@
 classdef Canonical < otp.ascherlineardae.AscherLinearDAEProblem
-    % The original problem defined by Uri Ascher in :cite:p:`Asc89` with 
-    % $\beta = 0.5$
+    % The problem defined by Uri Ascher in :cite:p:`Asc89` (sec. 2) 
+    % which uses timespan $t \in [0, 1]$  and  $\beta = 1 $.
     % 
     methods
-        function obj = Canonical(beta)
-            tspan = [0.0; 1];
+        function obj = Canonical(varargin)
+            % Create the Canonical CUSP problem object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. The accepted names are
+            %
+            %    - ``Beta`` – The parameter of the Ascher linear DAE problem.
             
-            if nargin < 1
-                beta = 0.5;
-            end
+            p = inputParser;
+            p.addParameter('beta', 1);
+            p.parse(varargin{:});
+            opts = p.Results;
             
-            params = otp.ascherlineardae.AscherLinearDAEParameters;
-            params.Beta = beta;
+            params      = otp.ascherlineardae.AscherLinearDAEParameters;
+            params.Beta = opts.beta;
             
-            y0 = [1; beta];
+            y0    = [1; params.Beta];
+            tspan = [0.0; 1.0];
             
             obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
new file mode 100644
index 00000000..46903e43
--- /dev/null
+++ b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
@@ -0,0 +1,19 @@
+classdef Petzold < otp.ascherlineardae.AscherLinearDAEProblem
+    % The Petzold DAE example :cite:p:`Pet86` as a variant of the 
+    % Ascher linear DAE problem
+    % which uses timespan $t \in [0, 1]$  and  $\beta = 0 $.
+    % 
+    methods
+        function obj = Petzold
+            % Create the Petzold example of the Ascher linear DAE problem object.
+            
+            params      = otp.ascherlineardae.AscherLinearDAEParameters;
+            params.Beta = 0;
+            
+            y0    = [1; params.Beta];
+            tspan = [0.0; 1.0];
+            
+            obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
+        end
+    end
+end
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
new file mode 100644
index 00000000..9d0ee1f0
--- /dev/null
+++ b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
@@ -0,0 +1,19 @@
+classdef Stiff < otp.ascherlineardae.AscherLinearDAEProblem
+    % The Stiff example from :cite:p:`Asc89` (sec. 2. A variant of the 
+    % Ascher linear DAE problem
+    % which uses timespan $t \in [0, 1]$  and  $\beta = 100 $.
+    % 
+    methods
+        function obj = Stiff
+            % Create the stiff example of the Ascher linear DAE problem object.
+            
+            params      = otp.ascherlineardae.AscherLinearDAEParameters;
+            params.Beta = 100;
+            
+            y0    = [1; params.Beta];
+            tspan = [0.0; 1.0];
+            
+            obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
+        end
+    end
+end
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
index b50416c0..595bcf3b 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
@@ -1,8 +1,8 @@
 classdef AscherLinearDAEParameters
     % Parameters for Ascher Linear DAE problem 
     properties
-        % Beta is a scalar parameter in the linear model. It also affects the initial condition
-        % of the algebraic variable.
-        Beta %MATLAB ONLY: (1,1) {mustBeNumeric} = 0.5
+        % Beta is a scalar parameter in the linear model. It affects the stifness 
+        % of the problem.
+        Beta %MATLAB ONLY: (1,1) {mustBeNumeric} = 1
     end
 end
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
index 195a709c..003d79a3 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
@@ -2,35 +2,33 @@
     % A simple linear differential algebraic problem.
     %
     % The Ascher linear DAE Problem :cite:p:`Asc89` is an index-1 differential agebraic 
-    % equation given by $M(t) y' = A(t) y + q(t) $ where
+    % equation given by
     %
     % $$
-    % M=\left(\begin{array}{cc}
+    % \begin{bmatrix}
     % 1 & -t \\
     % 0 & 0
-    % \end{array}\right), \quad A=\left(\begin{array}{cc}
+    % \end{bmatrix} \begin{bmatrix} y'(t) \\ z'(t) \end{bmatrix} = \left[ \begin{array}{cc}
     % -1 & 1+t \\
     % \beta & -1-\beta t
-    % \end{array}\right), \quad {q}=\left(\begin{array}{c}
+    % \end{array}\right] \begin{bmatrix} y(t) \\ z(t) \end{bmatrix} + \begin{bmatrix}
     % 0 \\
     % \sin t
-    % \end{array}\right), 
+    % \end{bmatrix}.
     % $$
     %
-    % defined on timespan $t \in [0,1]$, and initial condition $y_0 = [1, \beta]^T$. The exact solution
-    % is given by 
+    %  The problem has the following closed-form solution when the initial condition $y(0) = 1 , z(0) = \beta$ is used:
     %
     % $$
-    % y = \begin{pmatrix}
+    % \begin{bmatrix} y(t)\\ z(t) \end{bmatrix} = \begin{bmatrix}
     % t \sin(t) + (1 + \beta  t) e^{-t}\\
     % \beta  e^{-t} + \sin(t)
-    % \end{pmatrix}.
+    % \end{bmatrix}.
     % $$
-    %
-    % Due to its stiffness and time-dependant mass
-    % matrix, this simple DAE problem can 
-    % become challenging to solve. This problem is introduced in :cite:p:`Asc89` 
-    % to study the convergence of implcit solvers applied to DAEs.
+    % This DAE problem 
+    % can be used to investigate
+    % the convergence of implcit time-stepping methods due to its stiffness and time-dependant mass
+    % matrix.
     %
     % Notes
     % -----
@@ -39,14 +37,15 @@
     % +---------------------+-----------------------------------------+
     % | Number of Variables | 2                                       |
     % +---------------------+-----------------------------------------+
-    % | Stiff               | typically, depending on $\beta$         |
+    % | Stiff               | possibly, depending on $\beta$          |
     % +---------------------+-----------------------------------------+
     %
     % Example
     % -------
-    % >>> problem = otp.ascherlineardae.presets.Canonical(0.1);
-    % >>> sol = problem.solve('MaxStep',1e-5);
-    % >>> problem.plotPhaseSpace(sol);
+    % >>> problem = otp.ascherlineardae.presets.Canonical('Beta', 50);
+    % >>> t = linspace(0,1,100);
+    % >>> sol = problem.solveExcatly(t);
+    % >>> problem.plot(sol);
 
     properties (Access = private, Constant)
         NumComps = 2
@@ -65,11 +64,6 @@
             %    The initial conditions.
             % parameters : AscherLinearDAEParameters
             %    The parameters.
-            %
-            % Returns
-            % -------
-            % obj : AscherLinearDAEProblem
-            %    The constructed problem.
             obj@otp.Problem('Ascher Linear DAE', 2, timeSpan, y0, parameters);
         end
     end
@@ -85,7 +79,8 @@ function onSettingsChanged(obj)
                 'MassSingular', 'yes');
         end
         
-        function y = internalSolveExactly(obj, t)
+        function sol = internalSolveExactly(obj, t)
+            sol = [];
             beta = obj.Parameters.Beta;
             if ~isequal(obj.Y0, [1; beta])
                 error('OTP:noExactSolution', ...
@@ -94,6 +89,8 @@ function onSettingsChanged(obj)
             
             y = [t .* sin(t) + (1 + beta * t) .* exp(-t); ...
                 beta * exp(-t) + sin(t)];
+            sol.x = t;
+            sol.y = y;
         end
         
         function sol = internalSolve(obj, varargin)

From cce57cbda080000afea4e73dcb2f7b1f0654ee48 Mon Sep 17 00:00:00 2001
From: Andrey A Popov <apopov@vt.edu>
Date: Tue, 19 Sep 2023 16:00:03 -0500
Subject: [PATCH 05/20] fixes based on past conversations

---
 toolbox/+otp/+lorenz96/+presets/Canonical.m  | 16 +++++-----
 toolbox/+otp/+lorenz96/+presets/PopovSandu.m | 31 ++++++++++----------
 toolbox/+otp/+lorenz96/Lorenz96Problem.m     | 25 +++++++---------
 3 files changed, 34 insertions(+), 38 deletions(-)

diff --git a/toolbox/+otp/+lorenz96/+presets/Canonical.m b/toolbox/+otp/+lorenz96/+presets/Canonical.m
index ac4c3f56..da71ab3b 100644
--- a/toolbox/+otp/+lorenz96/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz96/+presets/Canonical.m
@@ -1,7 +1,7 @@
 classdef Canonical < otp.lorenz96.Lorenz96Problem
     % Original Lorenz '96 preset presented in :cite:p:`Lor96`
     % which uses time span $t \in [0, 720]$, $N = 40$, $F=8$, and initial
-    % conditions of $y_i = 8$ for all $i$ except for $y_20=8.008$. 
+    % conditions of $y_i = 8$ for all $i$ except for $y_{20}=8.008$. 
 
     methods
         function obj = Canonical(varargin)
@@ -9,15 +9,13 @@
             %
             % Parameters
             % ----------
-            % Size : numeric(1, 1)
-            %    The size of the problem as a positive integer.
-            % Forcing : numeric(N, 1)
-            %    The forcing as a vector of N constants.
+            % varargin
+            %    A variable number of name-value pairs. The accepted names are
+            %
+            %    - ``Size`` – The size of the problem as a positive integer.
+            %    - ``Forcing`` – The forcing as a scalar, vector of N constants, or as a
+            %    function.
             %
-            % Returns
-            % -------
-            % obj : Canonial
-            %    The constructed problem.
 
             p = inputParser;
             p.addParameter('Size', 40, @isscalar);
diff --git a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
index 3c171f1d..9f829934 100644
--- a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
+++ b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
@@ -1,9 +1,16 @@
 classdef PopovSandu < otp.lorenz96.Lorenz96Problem
     % A preset that has a cyclic forcing function that is different for
     % every variable. This preset was created for :cite:p:`PS19`.
-    % This preset uses time span $t \in [0, 720]$, $N = 40$, $F=8$, 
-    % four partitions, a forcing period of one time unit, and initial
-    % conditions of $y_i = 8$ for all $i$ except for $y_20=8.008$. 
+    % This preset uses time span $t \in [0, 720]$, $N = 40$, and initial
+    % conditions of $y_i = 8$ for all $i$ except for $y_{20}=8.008$. 
+    % The forcing, as a function of time is given by
+    %
+    % $$
+    % f(t) = 8 + 4\cos(2 \pi \omega (t+\text{mod}(i - 1, q)/q))
+    % $$
+    %
+    % where by default $q=4$ is the number of partitions, and $\omega = 1$
+    % is a forcing period of one time unit.
     
     methods
         function obj = PopovSandu(varargin)
@@ -11,19 +18,13 @@
             %
             % Parameters
             % ----------
-            % Size : numeric(1, 1)
-            %    The size of the problem as a positive integer.
-            % Partitions : numeric(1, 1)
-            %    The number of partitions into which to divide the
-            %    variables.
-            % ForcingPeriod : numeric(1, 1)
-            %    The period of the forcing function in radians per unit
-            %    time.
+            % varargin
+            %    A variable number of name-value pairs. The accepted names are
+            %
+            %    - ``Size`` – The size of the problem as a positive integer.
+            %    - ``Partitions`` – The number of partitions into which to divide the variables.
+            %    - ``ForcingPeriod`` – The period of the forcing function in radians per unit time.
             %
-            % Returns
-            % -------
-            % obj : PopovSandu
-            %    The constructed problem.
             
             p = inputParser;
             
diff --git a/toolbox/+otp/+lorenz96/Lorenz96Problem.m b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
index 72b3fd93..6bc62d6f 100644
--- a/toolbox/+otp/+lorenz96/Lorenz96Problem.m
+++ b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
@@ -1,22 +1,23 @@
 classdef Lorenz96Problem <  otp.Problem
-    % A chaotic system modeling the transfer of some quantity along some
-    % longitude.
-    %
-    %
+    % A chaotic system modeling nonlinear transfer of a dimensionless
+    % qauntity along a cyclic one dimensional domain.
+    % 
     % The $N$ variable dynamics :cite:p:`Lor96` are represented by the equation,
     %
     % $$
-    % y_i' = -y_{i-1}\left(y_{i-2} - y_{i+1}\right) - y_i + f(t),\quad i \in \mathbb{Z}_N,
+    % y_i' = -y_{i-1}\left(y_{i-2} - y_{i+1}\right) - y_i + f(t), \quad i = 1, \dots, N,
     % $$
-    %
-    % that exhibits chaotic behavior for certain values of the forcing function $f$.
+    % 
+    % where $y_0 = y_N$, $y_{-1} = y_{N - 1}$, and $y_{N + 1} = y_2$, exhibits 
+    % chaotic behavior for certain pairs of values of the dimension $N$ and
+    % forcing function $f$.
     %
     % Notes
     % -----
     % +---------------------+-----------------------------------------------------------+
     % | Type                | ODE                                                       |
     % +---------------------+-----------------------------------------------------------+
-    % | Number of Variables | $N$                                                       |
+    % | Number of Variables | $N$ for any postive integer four or greater               |
     % +---------------------+-----------------------------------------------------------+
     % | Stiff               | no                                                        |
     % +---------------------+-----------------------------------------------------------+
@@ -25,7 +26,7 @@
     % -------
     % >>> problem = otp.lorenz96.presets.Canonical('Forcing', @(t) 8 + 4*sin(t));
     % >>> sol = problem.solve();
-    % >>> problem.plotPhaseSpace(sol);
+    % >>> problem.movie(sol);
     %
     % See also
     % --------
@@ -43,15 +44,11 @@
             % ----------
             % timeSpan : numeric(1, 2)
             %    The start and final time.
-            % y0 : numeric(N, 1)
+            % y0 : numeric(:, 1)
             %    The initial conditions.
             % parameters : Lorenz96Parameters
             %    The parameters.
             %
-            % Returns
-            % -------
-            % obj : Lorenz96Problem
-            %    The constructed problem.
             obj@otp.Problem('Lorenz 96', [], timeSpan, y0, parameters);
         end
     end

From 0a65855d2eba99f2ca55af4a67966f971707b966 Mon Sep 17 00:00:00 2001
From: Andrey A Popov <apopov@vt.edu>
Date: Tue, 19 Sep 2023 20:37:19 -0500
Subject: [PATCH 06/20] constructor for lorenz 63 presets

---
 toolbox/+otp/+lorenz63/+presets/Canonical.m  | 4 ----
 toolbox/+otp/+lorenz63/+presets/LimitCycle.m | 3 +++
 toolbox/+otp/+lorenz63/+presets/Surprise.m   | 2 ++
 3 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/toolbox/+otp/+lorenz63/+presets/Canonical.m b/toolbox/+otp/+lorenz63/+presets/Canonical.m
index 74ba5ba0..bbdcae9e 100644
--- a/toolbox/+otp/+lorenz63/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz63/+presets/Canonical.m
@@ -16,10 +16,6 @@
             %    - ``rho`` – Value of $\rho$.
             %    - ``beta`` – Value of $\beta$.
             %
-            % Returns
-            % -------
-            % obj : Lorenz63Problem
-            %    The constructed problem.
 
             p = inputParser;
             p.addParameter('sigma', 10);
diff --git a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
index 1453b675..64076d11 100644
--- a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
+++ b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
@@ -5,6 +5,9 @@
     
     methods
         function obj = LimitCycle
+            % Create the LimitCycle Lorenz '63 problem object.
+            %
+            
             sigma = 10;
             rho   = 350;
             beta  = 8/3;
diff --git a/toolbox/+otp/+lorenz63/+presets/Surprise.m b/toolbox/+otp/+lorenz63/+presets/Surprise.m
index 27799f85..f0844cf2 100644
--- a/toolbox/+otp/+lorenz63/+presets/Surprise.m
+++ b/toolbox/+otp/+lorenz63/+presets/Surprise.m
@@ -4,6 +4,8 @@
 
     methods
         function obj = Surprise
+            % Create the Surprise Lorenz '63 problem object.
+            %
 
             sigma = 10;
             rho   = 100;

From 232afd0abb936f51b330659c0fafab250dad5b3e Mon Sep 17 00:00:00 2001
From: Andrey A Popov <apopov@vt.edu>
Date: Tue, 19 Sep 2023 21:05:26 -0500
Subject: [PATCH 07/20] fixes

---
 toolbox/+otp/+lorenz96/+presets/Canonical.m  | 2 +-
 toolbox/+otp/+lorenz96/+presets/PopovSandu.m | 2 +-
 toolbox/+otp/+lorenz96/Lorenz96Problem.m     | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/toolbox/+otp/+lorenz96/+presets/Canonical.m b/toolbox/+otp/+lorenz96/+presets/Canonical.m
index da71ab3b..c4ee68f8 100644
--- a/toolbox/+otp/+lorenz96/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz96/+presets/Canonical.m
@@ -1,7 +1,7 @@
 classdef Canonical < otp.lorenz96.Lorenz96Problem
     % Original Lorenz '96 preset presented in :cite:p:`Lor96`
     % which uses time span $t \in [0, 720]$, $N = 40$, $F=8$, and initial
-    % conditions of $y_i = 8$ for all $i$ except for $y_{20}=8.008$. 
+    % conditions of $y_i = 8$ for all $i$ except for $y_{\lfloor N/2 \rfloor}=8.008$. 
 
     methods
         function obj = Canonical(varargin)
diff --git a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
index 9f829934..2b7b4f4f 100644
--- a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
+++ b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
@@ -2,7 +2,7 @@
     % A preset that has a cyclic forcing function that is different for
     % every variable. This preset was created for :cite:p:`PS19`.
     % This preset uses time span $t \in [0, 720]$, $N = 40$, and initial
-    % conditions of $y_i = 8$ for all $i$ except for $y_{20}=8.008$. 
+    % conditions of $y_i = 8$ for all $i$ except for $y_{\lfloor N/2 \rfloor}=8.008$. 
     % The forcing, as a function of time is given by
     %
     % $$
diff --git a/toolbox/+otp/+lorenz96/Lorenz96Problem.m b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
index 6bc62d6f..459e9de8 100644
--- a/toolbox/+otp/+lorenz96/Lorenz96Problem.m
+++ b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
@@ -1,6 +1,6 @@
 classdef Lorenz96Problem <  otp.Problem
     % A chaotic system modeling nonlinear transfer of a dimensionless
-    % qauntity along a cyclic one dimensional domain.
+    % quantity along a cyclic one dimensional domain.
     % 
     % The $N$ variable dynamics :cite:p:`Lor96` are represented by the equation,
     %

From 5923c8b3bc35511878b17c510e9a02f49673c708 Mon Sep 17 00:00:00 2001
From: Andrey A Popov <apopov@vt.edu>
Date: Tue, 19 Sep 2023 21:20:42 -0500
Subject: [PATCH 08/20] changed variable names to unicode

---
 toolbox/+otp/+lorenz63/+presets/Canonical.m  | 10 +++++-----
 toolbox/+otp/+lorenz63/+presets/LimitCycle.m |  6 +++---
 toolbox/+otp/+lorenz63/+presets/Surprise.m   |  2 +-
 toolbox/+otp/+lorenz63/Lorenz63Problem.m     |  8 ++++----
 4 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/toolbox/+otp/+lorenz63/+presets/Canonical.m b/toolbox/+otp/+lorenz63/+presets/Canonical.m
index bbdcae9e..4a8b36e4 100644
--- a/toolbox/+otp/+lorenz63/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz63/+presets/Canonical.m
@@ -1,7 +1,7 @@
 classdef Canonical < otp.lorenz63.Lorenz63Problem
     % Original Lorenz '63 preset presented in :cite:p:`Lor63`
-    % which uses time span $t \in [0, 60]$, $\sigma = 10$, $\rho = 28$, 
-    % $\beta = 8/3$, and intial conditions $y_0 = [0, 1, 0]^T$.
+    % which uses time span $t \in [0, 60]$, $σ = 10$, $ρ = 28$, 
+    % $β = 8/3$, and intial conditions $y_0 = [0, 1, 0]^T$.
 
     methods
         function obj = Canonical(varargin)
@@ -12,9 +12,9 @@
             % varargin
             %    A variable number of name-value pairs. The accepted names are
             %
-            %    - ``sigma`` – Value of $\sigma$.
-            %    - ``rho`` – Value of $\rho$.
-            %    - ``beta`` – Value of $\beta$.
+            %    - ``sigma`` – Value of $σ$.
+            %    - ``rho`` – Value of $ρ$.
+            %    - ``beta`` – Value of $β$.
             %
 
             p = inputParser;
diff --git a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
index 64076d11..3f525411 100644
--- a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
+++ b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
@@ -1,13 +1,13 @@
 classdef LimitCycle < otp.lorenz63.Lorenz63Problem
     % Lorenz '63 preset limit cycle, a non-chaotic preset, from :cite:p:`Str18` 
-    % which uses time span $t \in [0, 60]$, $\sigma = 10$, $\rho = 350$, 
-    % $\beta = 8/3$, and intial conditions $y_0 = [0, 1, 0]^T$.
+    % which uses time span $t \in [0, 60]$, $σ = 10$, $ρ = 350$, 
+    % $β = 8/3$, and intial conditions $y_0 = [0, 1, 0]^T$.
     
     methods
         function obj = LimitCycle
             % Create the LimitCycle Lorenz '63 problem object.
             %
-            
+
             sigma = 10;
             rho   = 350;
             beta  = 8/3;
diff --git a/toolbox/+otp/+lorenz63/+presets/Surprise.m b/toolbox/+otp/+lorenz63/+presets/Surprise.m
index f0844cf2..58f383ad 100644
--- a/toolbox/+otp/+lorenz63/+presets/Surprise.m
+++ b/toolbox/+otp/+lorenz63/+presets/Surprise.m
@@ -1,6 +1,6 @@
 classdef Surprise < otp.lorenz63.Lorenz63Problem
     % Lorenz '63 preset 'surprise' from :cite:p:`Str18` which uses time span $t \in [0, 60]$, 
-    % $\sigma = 10$, $\rho = 100$, $\beta = 8/3$, and intial conditions $y_0 = [2, 1, 1]^T$.
+    % $σ = 10$, $ρ = 100$, $β = 8/3$, and intial conditions $y_0 = [2, 1, 1]^T$.
 
     methods
         function obj = Surprise
diff --git a/toolbox/+otp/+lorenz63/Lorenz63Problem.m b/toolbox/+otp/+lorenz63/Lorenz63Problem.m
index a60dace1..a9c41122 100644
--- a/toolbox/+otp/+lorenz63/Lorenz63Problem.m
+++ b/toolbox/+otp/+lorenz63/Lorenz63Problem.m
@@ -4,9 +4,9 @@
     % The three variable Lorenz '63 problem :cite:p:`Lor63` of the form,
     %
     % $$
-    % x' &= \sigma(y - x),\\
-    % y' &= \rho x - y - xz,\\
-    % z' &= xy - \beta z,
+    % x' &= σ(y - x),\\
+    % y' &= ρx - y - xz,\\
+    % z' &= xy - βz,
     % $$
     %
     % exhibits chaotic behavior for certain values of the parameters.
@@ -23,7 +23,7 @@
     % +---------------------+-----------------------------------------------------------+
     % | Number of Variables | 3                                                         |
     % +---------------------+-----------------------------------------------------------+
-    % | Stiff               | not typically, depending on $\sigma$, $\rho$, and $\beta$ |
+    % | Stiff               | not typically, depending on $σ$, $ρ$, and $β$             |
     % +---------------------+-----------------------------------------------------------+
     %
     % Example

From e4e4ebb44453b4c435213a26e77c31225fb7b7d7 Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Tue, 26 Sep 2023 21:10:10 -0700
Subject: [PATCH 09/20] Sort presents alphabetically

---
 docs/generate_problem_rst.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/generate_problem_rst.py b/docs/generate_problem_rst.py
index 1715723e..ae455dec 100644
--- a/docs/generate_problem_rst.py
+++ b/docs/generate_problem_rst.py
@@ -30,7 +30,7 @@
 .. automodule:: +otp.{problem.parent.name}.+presets
 ''')
 
-        for preset in problem.parent.glob('+presets/*.m'):
+        for preset in sorted(problem.parent.glob('+presets/*.m')):
             stream.write(f'''
 .. autoclass:: {preset.stem}
     :members:

From be1079b7787b8fa8c9d04f5797b89ffa6259c5b9 Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Sun, 1 Oct 2023 21:17:02 -0700
Subject: [PATCH 10/20] Feature/oregonator docs (#119)

* Start Oregonator docs

* Add Canonical constructor docs

* Fix Jacobian
---
 docs/references.bib                           | 35 ++++++++
 toolbox/+otp/+oregonator/+presets/Canonical.m | 41 +++++++++-
 .../+otp/+oregonator/+presets/EnrightHull.m   | 35 ++++++++
 .../+oregonator/+presets/GottwaldWanner.m     | 37 +++++++++
 .../+otp/+oregonator/OregonatorParameters.m   | 15 +++-
 toolbox/+otp/+oregonator/OregonatorProblem.m  | 79 +++++++++++++++++--
 toolbox/+otp/+oregonator/f.m                  |  8 +-
 toolbox/+otp/+oregonator/jacobian.m           |  8 +-
 8 files changed, 242 insertions(+), 16 deletions(-)
 create mode 100644 toolbox/+otp/+oregonator/+presets/EnrightHull.m
 create mode 100644 toolbox/+otp/+oregonator/+presets/GottwaldWanner.m

diff --git a/docs/references.bib b/docs/references.bib
index af3d0a2f..be483b64 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -120,3 +120,38 @@ @article{PS19
   publisher={Copernicus Publications G{\"o}ttingen, Germany},
   doi={10.5194/npg-26-109-2019}
 }
+
+@article{FN74,
+    author = {Field, Richard J. and Noyes, Richard M.},
+    title = "{Oscillations in chemical systems. {IV}. {L}imit cycle behavior in a model of a real chemical reaction}",
+    journal = {The Journal of Chemical Physics},
+    volume = {60},
+    number = {5},
+    pages = {1877-1884},
+    year = {1974},
+    issn = {0021-9606},
+    doi = {10.1063/1.1681288}
+}
+
+@incollection{EH76,
+    title = {Comparing Numerical Methods for the Solution of Stiff Systems of {ODE}s Arising in Chemistry},
+    editor = {L. LAPIDUS and W.E. SCHIESSER},
+    booktitle = {Numerical Methods for Differential Systems},
+    publisher = {Academic Press},
+    pages = {45-66},
+    year = {1976},
+    isbn = {978-0-12-436640-4},
+    doi = {10.1016/B978-0-12-436640-4.50008-3},
+    author = {W.H. Enright and T.E. Hull}
+}
+
+@article{GW82,
+    author = {Björn A. Gottwald and Gerhard Wanner},
+    title ={Comparison of Numerical Methods for Stiff Differential Equations in Biology and Chemistry},
+    journal = {SIMULATION},
+    volume = {38},
+    number = {2},
+    pages = {61-66},
+    year = {1982},
+    doi = {10.1177/003754978203800206}
+}
diff --git a/toolbox/+otp/+oregonator/+presets/Canonical.m b/toolbox/+otp/+oregonator/+presets/Canonical.m
index fa6bd123..7c40b277 100644
--- a/toolbox/+otp/+oregonator/+presets/Canonical.m
+++ b/toolbox/+otp/+oregonator/+presets/Canonical.m
@@ -1,9 +1,48 @@
 classdef Canonical < otp.oregonator.OregonatorProblem
+    % The Oregonator configuration used in :cite:p:`HW96` (p. 144) called OREGO. It uses time span $t \in [0, 360]$,
+    % initial condition $y_0 = [1, 2, 3]^T$, and parameters
+    %
+    % $$
+    % f &= 1, \\
+    % q &= 8.375 \cdot 10^{-6}, \\
+    % s &= 77.27, \\
+    % w &= 0.1610.
+    % $$
+
     methods
-        function obj = Canonical
+        function obj = Canonical(varargin)
+            % Create the Canonical problem object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. The accepted names are
+            %
+            %    - ``f`` – Value of $f$.
+            %    - ``q`` – Value of $q$.
+            %    - ``s`` – Value of $s$.
+            %    - ``w`` – Value of $w$.
+            %
+            % Returns
+            % -------
+            % obj : Canonical
+            %    The constructed problem.
+
+            p = inputParser;
+            p.addParameter('f', 1);
+            p.addParameter('q', 8.375e-6);
+            p.addParameter('s', 77.27);
+            p.addParameter('w', 0.1610);
+            p.parse(varargin{:});
+            opts = p.Results;
+
             tspan = [0, 360];
             y0 = [1; 2; 3];
             params = otp.oregonator.OregonatorParameters;
+            params.F = 1;
+            params.Q = 8.375e-6;
+            params.S = 77.27;
+            params.W = 0.1610;
             
             obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+oregonator/+presets/EnrightHull.m b/toolbox/+otp/+oregonator/+presets/EnrightHull.m
new file mode 100644
index 00000000..38eff49b
--- /dev/null
+++ b/toolbox/+otp/+oregonator/+presets/EnrightHull.m
@@ -0,0 +1,35 @@
+classdef EnrightHull < otp.oregonator.OregonatorProblem
+    % The Oregonator configuration used in problem CHM 9 of :cite:p:`EH76`. It uses time span $t \in [0, 300]$, initial
+    % condition $y_0 = [4, 1.1, 4]^T$, and parameters
+    %
+    % $$
+    % f &= 1, \\
+    % q &= 8.375 \cdot 10^{-6}, \\
+    % s &= 77.27, \\
+    % w &= 0.1610.
+    % $$
+
+    methods
+        function obj = EnrightHull
+            % Create the EnrightHull problem object.
+            %
+            % Parameters
+            % ----------
+            %
+            % Returns
+            % -------
+            % obj : EnrightHull
+            %    The constructed problem.
+
+            tspan = [0, 300];
+            y0 = [4; 1.1; 4];
+            params = otp.oregonator.OregonatorParameters;
+            params.F = 1;
+            params.Q = 8.375e-6;
+            params.S = 77.27;
+            params.W = 0.1610;
+            
+            obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
+        end
+    end
+end
\ No newline at end of file
diff --git a/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m b/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
new file mode 100644
index 00000000..3dfc4c9b
--- /dev/null
+++ b/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
@@ -0,0 +1,37 @@
+classdef GottwaldWanner < otp.oregonator.OregonatorProblem
+    % The Oregonator configuration from :cite:p:`GW82` which uses time span $t \in [0, 302.85805]$, initial condition
+    % $y_0 = [4, 1.331391, 2.852348]^T$, and parameters
+    %
+    % $$
+    % f &= 1, \\
+    % q &= 8.375 \cdot 10^{-6}, \\
+    % s &= 77.27, \\
+    % w &= 0.1610.
+    % $$
+    %
+    % The initial condition lies on a stable limit cycle, and the time span is chosen to be one period.
+
+    methods
+        function obj = GottwaldWanner
+            % Create the GottwaldWanner problem object.
+            %
+            % Parameters
+            % ----------
+            %
+            % Returns
+            % -------
+            % obj : GottwaldWanner
+            %    The constructed problem.
+
+            tspan = [0, 302.85805];
+            y0 = [4; 1.331391; 2.852348];
+            params = otp.oregonator.OregonatorParameters;
+            params.F = 1;
+            params.Q = 8.375e-6;
+            params.S = 77.27;
+            params.W = 0.1610;
+            
+            obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
+        end
+    end
+end
\ No newline at end of file
diff --git a/toolbox/+otp/+oregonator/OregonatorParameters.m b/toolbox/+otp/+oregonator/OregonatorParameters.m
index 87ba9d5d..ff354805 100644
--- a/toolbox/+otp/+oregonator/OregonatorParameters.m
+++ b/toolbox/+otp/+oregonator/OregonatorParameters.m
@@ -1,4 +1,17 @@
 classdef OregonatorParameters
-    %OREGONATORPARAMETERS
+    % Parameters for the Oregonator problem.
+    properties
+        % The stoichiometric factor $f$.
+        F %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
+        
+        % The reaction constant $q$.
+        Q %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
+        
+        % The reaction constant $s$.
+        S %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
+        
+        % The reaction constant $w$.
+        W %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
+    end
 end
 
diff --git a/toolbox/+otp/+oregonator/OregonatorProblem.m b/toolbox/+otp/+oregonator/OregonatorProblem.m
index 7126f048..d4c57e10 100644
--- a/toolbox/+otp/+oregonator/OregonatorProblem.m
+++ b/toolbox/+otp/+oregonator/OregonatorProblem.m
@@ -1,25 +1,92 @@
 classdef OregonatorProblem < otp.Problem
+    % A periodic, three-variable model for the Belousov–Zhabotinsky reaction.
+    %
+    % The Oregonator chemical reaction :cite:p:`FN74` is given by
+    %
+    %
+    % $$
+    % \ce{
+    % A + Y &-> X + P \\
+    % X + Y &-> 2P \\
+    % A + X &-> 2X + 2Z \\
+    % 2X &-> A + P \\
+    % B + Z &-> 1/2 $f$ Y.
+    % }
+    % $$
+    %
+    % Species $\ce{A = BrO3-}$, $\ce{B = CH2(COOH)2}$, and $\ce{P = HOBr}$ or $\ce{BrCH(COOH)2}$ are assumed to be
+    % constant. The dynamic concentrations of intermediate species $\ce{X = HBrO2}$, $\ce{Y = Br-}$, and
+    % $\ce{Z = Ce(IV)}$ can be modeled by an ODE in three variables. Field and Noyes :cite:p:`FN74` proposed the
+    % nondimensionalized form
+    %
+    % $$
+    % α' &= s(η - η α + α - q α^2), \\
+    % η' &= s^{-1}(-η - η α + f ρ), \\
+    % ρ' &= w (α - ρ),
+    % $$
+    %
+    % where $α$, $η$, and $ρ$ are scaled concentrations of $\ce{X}$, $\ce{Y}$, and $\ce{Z}$, respectively.
+    %
+    % Notes
+    % -----
+    % +---------------------+------------------------------------------------+
+    % | Type                | ODE                                            |
+    % +---------------------+------------------------------------------------+
+    % | Number of Variables | 3                                              |
+    % +---------------------+------------------------------------------------+
+    % | Stiff               | typically, depending on $f$, $q$, $s$, and $w$ |
+    % +---------------------+------------------------------------------------+
+    %
+    % Example
+    % -------
+    % >>> problem = otp.oregonator.presets.Canonical;
+    % >>> sol = problem.solve();
+    % >>> problem.plotPhaseSpace(sol, 'Vars', [2, 3]);
+
     methods
         function obj = OregonatorProblem(timeSpan, y0, parameters)
+            % Create a Oregonator problem object.
+            %
+            % Parameters
+            % ----------
+            % timeSpan : numeric(1, 2)
+            %    The start and final time.
+            % y0 : numeric(3, 1)
+            %    The initial conditions.
+            % parameters : OregonatorParameters
+            %    The parameters.
+            %
+            % Returns
+            % -------
+            % obj : OregonatorProblem
+            %    The constructed problem.
             obj@otp.Problem('Oregonator', 3, timeSpan, y0, parameters);
         end
     end
     
     methods (Access = protected)
         function onSettingsChanged(obj)
-            obj.RHS = otp.RHS(@otp.oregonator.f, ...
-                'Jacobian', @otp.oregonator.jacobian, ...
+            f = obj.Parameters.F;
+            q = obj.Parameters.Q;
+            s = obj.Parameters.S;
+            w = obj.Parameters.W;
+            
+            obj.RHS = otp.RHS(@(t, y) otp.oregonator.f(t, y, f, q, s, w), ...
+                'Jacobian', @(t, y) otp.oregonator.jacobian(t, y, f, q, s, w), ...
                 'Vectorized', 'on');
         end
         
         function fig = internalPlot(obj, t, y, varargin)
-            fig = internalPlot@otp.Problem(obj, t, y, ...
-                'yscale', 'log', varargin{:});
+            fig = internalPlot@otp.Problem(obj, t, y, 'yscale', 'log', varargin{:});
+        end
+        
+        function fig = internalPlotPhaseSpace(obj, t, y, varargin)
+            fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, 'xscale', 'log', 'yscale', 'log', 'zscale', 'log', ...
+                varargin{:});
         end
         
         function mov = internalMovie(obj, t, y, varargin)
-            mov = internalMovie@otp.Problem(obj, t, y, ...
-                'yscale', 'log', varargin{:});
+            mov = internalMovie@otp.Problem(obj, t, y, 'yscale', 'log', varargin{:});
         end
     end
 end
diff --git a/toolbox/+otp/+oregonator/f.m b/toolbox/+otp/+oregonator/f.m
index 98b7a787..ebc53002 100644
--- a/toolbox/+otp/+oregonator/f.m
+++ b/toolbox/+otp/+oregonator/f.m
@@ -1,12 +1,12 @@
-function dy = f(~, y)
+function dy = f(~, y, f, q, s, w)
 
 y1 = y(1, :);
 y2 = y(2, :);
 y3 = y(3, :);
 
 dy = [ ...
-    77.27 * (y2 + y1 .* (1 - 8.375e-6 * y1 - y2)); ...
-    (y3 - (1 + y1) .* y2) / 77.27; ...
-    0.161 * (y1 - y3)];
+    s * (y2 + y1 .* (1 - q * y1 - y2)); ...
+    (f * y3 - (1 + y1) .* y2) / s; ...
+    w * (y1 - y3)];
 
 end
diff --git a/toolbox/+otp/+oregonator/jacobian.m b/toolbox/+otp/+oregonator/jacobian.m
index ea9f3fc4..ad3992df 100644
--- a/toolbox/+otp/+oregonator/jacobian.m
+++ b/toolbox/+otp/+oregonator/jacobian.m
@@ -1,11 +1,11 @@
-function J = jacobian(~, y)
+function J = jacobian(~, y, f, q, s, w)
 
 y1 = y(1);
 y2 = y(2);
 
 J = [ ...
-    77.27 * (1 - y2 - 1.675e-5 * y1), 77.27 * (1 - y1), 0; ...
-    -y2 / 77.27, -(1 + y1) / 77.27, 1 / 77.27; ...
-    0.161, 0, -0.161];
+    s * (1 - 2 * q * y1 - y2), s * (1 - y1), 0; ...
+    -y2 / s, -(1 + y1) / s, f / s; ...
+    w, 0, -w];
 
 end

From 1765fe31bd8a258cfbe411d69d35aff6bfaf203c Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Sun, 1 Oct 2023 21:57:43 -0700
Subject: [PATCH 11/20] Fix sphinx warning

---
 toolbox/+otp/+lorenz96/+presets/Canonical.m | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/toolbox/+otp/+lorenz96/+presets/Canonical.m b/toolbox/+otp/+lorenz96/+presets/Canonical.m
index c4ee68f8..a21c4784 100644
--- a/toolbox/+otp/+lorenz96/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz96/+presets/Canonical.m
@@ -13,8 +13,7 @@
             %    A variable number of name-value pairs. The accepted names are
             %
             %    - ``Size`` – The size of the problem as a positive integer.
-            %    - ``Forcing`` – The forcing as a scalar, vector of N constants, or as a
-            %    function.
+            %    - ``Forcing`` – The forcing as a scalar, vector of N constants, or as a function.
             %
 
             p = inputParser;

From fb2d7e71431499cc20a690ac16bd3b2874027885 Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Sun, 1 Oct 2023 22:56:41 -0700
Subject: [PATCH 12/20] Update python packages

---
 docs/requirements.txt | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/requirements.txt b/docs/requirements.txt
index 8f463ca3..26f58a86 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -1,9 +1,9 @@
 sphinx==6.2.1
 sphinx-autobuild==2021.3.14
-sphinx_rtd_theme==1.2.0
-sphinxcontrib-matlabdomain==0.18.0
+sphinx_rtd_theme==1.3.0
+sphinxcontrib-matlabdomain==0.20.2
 sphinxcontrib-napoleon==0.7
 sphinx-math-dollar==1.2.1
 sphinxcontrib-bibtex==2.5.0
-myst-parser==1.0.0
-pyyaml==6.0
+myst-parser==2.0.0
+pyyaml==6.0.1

From 3cfaec9e06fa438fb915787c9d50f1d5a8446dff Mon Sep 17 00:00:00 2001
From: Arash Sarshar <elswit@users.noreply.github.com>
Date: Tue, 3 Oct 2023 16:54:59 -0700
Subject: [PATCH 13/20] up Ascher docs re code review

---
 .../+ascherlineardae/+presets/Canonical.m     |  2 +-
 .../+otp/+ascherlineardae/+presets/Petzold.m  |  6 ++--
 .../+otp/+ascherlineardae/+presets/Stiff.m    |  5 ++-
 .../+ascherlineardae/AscherLinearDAEProblem.m | 34 +++++++++----------
 4 files changed, 23 insertions(+), 24 deletions(-)

diff --git a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
index 865a3c94..e47f7028 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
@@ -1,6 +1,6 @@
 classdef Canonical < otp.ascherlineardae.AscherLinearDAEProblem
     % The problem defined by Uri Ascher in :cite:p:`Asc89` (sec. 2) 
-    % which uses timespan $t \in [0, 1]$  and  $\beta = 1 $.
+    % which uses timespan $t \in [0, 1]$  and intial condition $[y_0, z_0]^T = [1, \beta]^T  $.
     % 
     methods
         function obj = Canonical(varargin)
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
index 46903e43..1aba156d 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
@@ -1,7 +1,7 @@
 classdef Petzold < otp.ascherlineardae.AscherLinearDAEProblem
-    % The Petzold DAE example :cite:p:`Pet86` as a variant of the 
-    % Ascher linear DAE problem
-    % which uses timespan $t \in [0, 1]$  and  $\beta = 0 $.
+    % The Petzold DAE example :cite:p:`Pet86` as a special case of the 
+    % Ascher linear DAE problem.
+    % This preset uses timespan $t \in [0, 1]$  and  $\beta = 0 $ with the initial condition $[y_0, z_0]^T = [1, 0]^T $.
     % 
     methods
         function obj = Petzold
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
index 9d0ee1f0..7cfdaac4 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
@@ -1,12 +1,11 @@
 classdef Stiff < otp.ascherlineardae.AscherLinearDAEProblem
-    % The Stiff example from :cite:p:`Asc89` (sec. 2. A variant of the 
+    % The Stiff example from :cite:p:`Asc89`. A variant of the 
     % Ascher linear DAE problem
-    % which uses timespan $t \in [0, 1]$  and  $\beta = 100 $.
+    % which uses timespan $t \in [0, 1]$  and  $\beta = 100 $ with the initial condition $[y_0, z_0]^T = [1, 100]^T $.
     % 
     methods
         function obj = Stiff
             % Create the stiff example of the Ascher linear DAE problem object.
-            
             params      = otp.ascherlineardae.AscherLinearDAEParameters;
             params.Beta = 100;
             
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
index 003d79a3..9cf20d82 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
@@ -1,7 +1,7 @@
 classdef AscherLinearDAEProblem < otp.Problem
-    % A simple linear differential algebraic problem.
+    % A linear differential-algebraic problem with time-dependant mass matrix.
     %
-    % The Ascher linear DAE Problem :cite:p:`Asc89` is an index-1 differential agebraic 
+    % The Ascher linear DAE Problem :cite:p:`Asc89` is an index-1 differential-agebraic 
     % equation given by
     %
     % $$
@@ -17,7 +17,7 @@
     % \end{bmatrix}.
     % $$
     %
-    %  The problem has the following closed-form solution when the initial condition $y(0) = 1 , z(0) = \beta$ is used:
+    % When the initial condition $y(0) = 1 , z(0) = \beta$ is used, the problem has the following closed-form solution:
     %
     % $$
     % \begin{bmatrix} y(t)\\ z(t) \end{bmatrix} = \begin{bmatrix}
@@ -43,14 +43,9 @@
     % Example
     % -------
     % >>> problem = otp.ascherlineardae.presets.Canonical('Beta', 50);
-    % >>> t = linspace(0,1,100);
-    % >>> sol = problem.solveExcatly(t);
-    % >>> problem.plot(sol);
-
-    properties (Access = private, Constant)
-        NumComps = 2
-        VarNames = 'yz'
-    end
+    % >>> t = linspace(0,1);
+    % >>> sol = problem.solveExactly(t);
+    % >>> problem.plot(t, sol);
 
     methods
         function obj = AscherLinearDAEProblem(timeSpan, y0, parameters)
@@ -61,7 +56,7 @@
             % timeSpan : numeric(1, 2)
             %    The start and final time.
             % y0 : numeric(2, 1)
-            %    The initial conditions.
+            %    The initial condition.
             % parameters : AscherLinearDAEParameters
             %    The parameters.
             obj@otp.Problem('Ascher Linear DAE', 2, timeSpan, y0, parameters);
@@ -78,9 +73,16 @@ function onSettingsChanged(obj)
                 'MStateDependence', 'none', ...
                 'MassSingular', 'yes');
         end
-        
-        function sol = internalSolveExactly(obj, t)
-            sol = [];
+
+        function label = internalIndex2label(~, index)
+            if index == 1
+                label = 'Differential Variable';
+            else
+                label = 'Algebraic Variable';
+            end
+        end
+
+        function y = internalSolveExactly(obj, t)
             beta = obj.Parameters.Beta;
             if ~isequal(obj.Y0, [1; beta])
                 error('OTP:noExactSolution', ...
@@ -89,8 +91,6 @@ function onSettingsChanged(obj)
             
             y = [t .* sin(t) + (1 + beta * t) .* exp(-t); ...
                 beta * exp(-t) + sin(t)];
-            sol.x = t;
-            sol.y = y;
         end
         
         function sol = internalSolve(obj, varargin)

From f05b7e7d849069f1e383f3099931b7f4d27cffde Mon Sep 17 00:00:00 2001
From: Arash Sarshar <elswit@users.noreply.github.com>
Date: Tue, 3 Oct 2023 19:55:28 -0700
Subject: [PATCH 14/20] up Ascher docs, line width and unicode characters fixed

---
 .../+ascherlineardae/+presets/Canonical.m     |  6 ++---
 .../+otp/+ascherlineardae/+presets/Petzold.m  |  5 ++--
 .../+otp/+ascherlineardae/+presets/Stiff.m    |  5 ++--
 .../AscherLinearDAEParameters.m               |  3 +--
 .../+ascherlineardae/AscherLinearDAEProblem.m | 27 +++++++++----------
 5 files changed, 20 insertions(+), 26 deletions(-)

diff --git a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
index e47f7028..c092d703 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
@@ -1,6 +1,6 @@
 classdef Canonical < otp.ascherlineardae.AscherLinearDAEProblem
-    % The problem defined by Uri Ascher in :cite:p:`Asc89` (sec. 2) 
-    % which uses timespan $t \in [0, 1]$  and intial condition $[y_0, z_0]^T = [1, \beta]^T  $.
+    % The problem defined by Uri Ascher in :cite:p:`Asc89` (sec. 2) which uses time span $t \in [0, 1]$ and intial 
+    % condition $[y_0, z_0]^T = [1, β]^T$.
     % 
     methods
         function obj = Canonical(varargin)
@@ -11,7 +11,7 @@
             % varargin
             %    A variable number of name-value pairs. The accepted names are
             %
-            %    - ``Beta`` – The parameter of the Ascher linear DAE problem.
+            %    - ``Beta`` – Value of $β$.
             
             p = inputParser;
             p.addParameter('beta', 1);
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
index 1aba156d..d8f4947d 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
@@ -1,7 +1,6 @@
 classdef Petzold < otp.ascherlineardae.AscherLinearDAEProblem
-    % The Petzold DAE example :cite:p:`Pet86` as a special case of the 
-    % Ascher linear DAE problem.
-    % This preset uses timespan $t \in [0, 1]$  and  $\beta = 0 $ with the initial condition $[y_0, z_0]^T = [1, 0]^T $.
+    % The Petzold DAE example :cite:p:`Pet86` as a special case of the Ascher linear DAE problem. This preset uses time
+    % span $t \in [0, 1]$ and $β = 0 $ with the initial condition $[y_0, z_0]^T = [1, 0]^T $.
     % 
     methods
         function obj = Petzold
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
index 7cfdaac4..c3ab60ce 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
@@ -1,7 +1,6 @@
 classdef Stiff < otp.ascherlineardae.AscherLinearDAEProblem
-    % The Stiff example from :cite:p:`Asc89`. A variant of the 
-    % Ascher linear DAE problem
-    % which uses timespan $t \in [0, 1]$  and  $\beta = 100 $ with the initial condition $[y_0, z_0]^T = [1, 100]^T $.
+    % The Stiff example from :cite:p:`Asc89`. A variant of the Ascher linear DAE problem which uses time span 
+    % $t \in [0, 1]$ and $β = 100$ with the initial condition $[y_0, z_0]^T = [1, 100]^T$.
     % 
     methods
         function obj = Stiff
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
index 595bcf3b..ca31dccc 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
@@ -1,8 +1,7 @@
 classdef AscherLinearDAEParameters
     % Parameters for Ascher Linear DAE problem 
     properties
-        % Beta is a scalar parameter in the linear model. It affects the stifness 
-        % of the problem.
+        % A scalar parameter $β$ in the linear model. It affects the stifness of the problem.
         Beta %MATLAB ONLY: (1,1) {mustBeNumeric} = 1
     end
 end
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
index 9cf20d82..e68c84af 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
@@ -1,8 +1,7 @@
 classdef AscherLinearDAEProblem < otp.Problem
     % A linear differential-algebraic problem with time-dependant mass matrix.
     %
-    % The Ascher linear DAE Problem :cite:p:`Asc89` is an index-1 differential-agebraic 
-    % equation given by
+    % The Ascher linear DAE Problem :cite:p:`Asc89` is an index-1 differential-agebraic equation given by
     %
     % $$
     % \begin{bmatrix}
@@ -10,25 +9,23 @@
     % 0 & 0
     % \end{bmatrix} \begin{bmatrix} y'(t) \\ z'(t) \end{bmatrix} = \left[ \begin{array}{cc}
     % -1 & 1+t \\
-    % \beta & -1-\beta t
+    % β & -1-β t
     % \end{array}\right] \begin{bmatrix} y(t) \\ z(t) \end{bmatrix} + \begin{bmatrix}
     % 0 \\
-    % \sin t
+    % \sin(t)
     % \end{bmatrix}.
     % $$
     %
-    % When the initial condition $y(0) = 1 , z(0) = \beta$ is used, the problem has the following closed-form solution:
+    % When the initial condition $y(0) = 1 , z(0) = β$ is used, the problem has the following closed-form solution:
     %
     % $$
     % \begin{bmatrix} y(t)\\ z(t) \end{bmatrix} = \begin{bmatrix}
-    % t \sin(t) + (1 + \beta  t) e^{-t}\\
-    % \beta  e^{-t} + \sin(t)
+    % t \sin(t) + (1 + β  t) e^{-t}\\
+    % β  e^{-t} + \sin(t)
     % \end{bmatrix}.
     % $$
-    % This DAE problem 
-    % can be used to investigate
-    % the convergence of implcit time-stepping methods due to its stiffness and time-dependant mass
-    % matrix.
+    % This DAE problem can be used to investigate the convergence of implcit time-stepping methods due to its stiffness
+    % and time-dependant mass matrix.
     %
     % Notes
     % -----
@@ -37,13 +34,13 @@
     % +---------------------+-----------------------------------------+
     % | Number of Variables | 2                                       |
     % +---------------------+-----------------------------------------+
-    % | Stiff               | possibly, depending on $\beta$          |
+    % | Stiff               | possibly, depending on $β$              |
     % +---------------------+-----------------------------------------+
     %
     % Example
     % -------
     % >>> problem = otp.ascherlineardae.presets.Canonical('Beta', 50);
-    % >>> t = linspace(0,1);
+    % >>> t = linspace(0, 1);
     % >>> sol = problem.solveExactly(t);
     % >>> problem.plot(t, sol);
 
@@ -76,9 +73,9 @@ function onSettingsChanged(obj)
 
         function label = internalIndex2label(~, index)
             if index == 1
-                label = 'Differential Variable';
+                label = 'Differential';
             else
-                label = 'Algebraic Variable';
+                label = 'Algebraic';
             end
         end
 

From 3e6ba2ae43c8b65ab7f798692d33fe091a5edf84 Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Tue, 17 Oct 2023 19:56:37 -0700
Subject: [PATCH 15/20] Make superclass for params (#118)

* Make superclass for params

* Fix param capitalization

* Add docs

* Update Ascher Linear DAE

* Update CUSP

* Update L'96

* Update Oregonator

* Update PR

* Update Robertson

* Typo fixes

* Fix Greek letters
---
 docs/structure/index.rst                      |  4 ++-
 docs/structure/parameters.rst                 |  6 ++++
 .../+ascherlineardae/+presets/Canonical.m     | 11 ++----
 .../+otp/+ascherlineardae/+presets/Petzold.m  |  6 ++--
 .../+otp/+ascherlineardae/+presets/Stiff.m    |  7 ++--
 .../AscherLinearDAEParameters.m               | 16 ++++++++-
 toolbox/+otp/+cusp/+presets/Canonical.m       | 25 ++++++-------
 toolbox/+otp/+cusp/CUSPParameters.m           | 20 +++++++++--
 toolbox/+otp/+cusp/CUSPProblem.m              | 31 ++++++++--------
 toolbox/+otp/+lorenz63/+presets/Canonical.m   | 15 +-------
 toolbox/+otp/+lorenz63/+presets/LimitCycle.m  | 11 +-----
 toolbox/+otp/+lorenz63/+presets/Surprise.m    | 12 +------
 toolbox/+otp/+lorenz63/Lorenz63Parameters.m   | 16 ++++++++-
 toolbox/+otp/+lorenz96/+presets/Canonical.m   | 19 ++++------
 toolbox/+otp/+lorenz96/Lorenz96Parameters.m   | 16 ++++++++-
 toolbox/+otp/+oregonator/+presets/Canonical.m | 14 +-------
 .../+otp/+oregonator/+presets/EnrightHull.m   |  6 +---
 .../+oregonator/+presets/GottwaldWanner.m     |  6 +---
 .../+otp/+oregonator/OregonatorParameters.m   | 16 ++++++++-
 .../+protherorobinson/+presets/Canonical.m    | 36 +++++++------------
 .../ProtheroRobinsonParameters.m              | 22 +++++++++---
 .../ProtheroRobinsonProblem.m                 | 18 +++++-----
 toolbox/+otp/+robertson/+presets/Canonical.m  | 14 ++++----
 toolbox/+otp/+robertson/+presets/ROBER.m      |  9 +----
 toolbox/+otp/+robertson/RobertsonParameters.m | 16 ++++++++-
 toolbox/+otp/Parameters.m                     | 28 +++++++++++++++
 toolbox/+otp/Problem.m                        |  2 +-
 toolbox/+otp/RHS.m                            |  2 +-
 28 files changed, 225 insertions(+), 179 deletions(-)
 create mode 100644 docs/structure/parameters.rst
 create mode 100644 toolbox/+otp/Parameters.m

diff --git a/docs/structure/index.rst b/docs/structure/index.rst
index b5798d8c..0e555eaa 100644
--- a/docs/structure/index.rst
+++ b/docs/structure/index.rst
@@ -8,4 +8,6 @@ TODO
    :maxdepth: 1
    :glob:
    
-   *
+   problem
+   rhs
+   parameters
diff --git a/docs/structure/parameters.rst b/docs/structure/parameters.rst
new file mode 100644
index 00000000..1cc3d570
--- /dev/null
+++ b/docs/structure/parameters.rst
@@ -0,0 +1,6 @@
+Parameters
+================================================================================
+
+.. automodule:: +otp
+.. autoclass:: Parameters
+    :members:
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
index c092d703..1f6a35c5 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
@@ -13,15 +13,8 @@
             %
             %    - ``Beta`` – Value of $β$.
             
-            p = inputParser;
-            p.addParameter('beta', 1);
-            p.parse(varargin{:});
-            opts = p.Results;
-            
-            params      = otp.ascherlineardae.AscherLinearDAEParameters;
-            params.Beta = opts.beta;
-            
-            y0    = [1; params.Beta];
+            params = otp.ascherlineardae.AscherLinearDAEParameters('Beta', 1, varargin{:});
+            y0 = [1; params.Beta];
             tspan = [0.0; 1.0];
             
             obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
index d8f4947d..0e67717e 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
@@ -6,10 +6,8 @@
         function obj = Petzold
             % Create the Petzold example of the Ascher linear DAE problem object.
             
-            params      = otp.ascherlineardae.AscherLinearDAEParameters;
-            params.Beta = 0;
-            
-            y0    = [1; params.Beta];
+            params = otp.ascherlineardae.AscherLinearDAEParameters('Beta', 0);
+            y0 = [1; params.Beta];
             tspan = [0.0; 1.0];
             
             obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
index c3ab60ce..a32cbfe7 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
@@ -5,10 +5,9 @@
     methods
         function obj = Stiff
             % Create the stiff example of the Ascher linear DAE problem object.
-            params      = otp.ascherlineardae.AscherLinearDAEParameters;
-            params.Beta = 100;
-            
-            y0    = [1; params.Beta];
+
+            params = otp.ascherlineardae.AscherLinearDAEParameters('Beta', 100);
+            y0 = [1; params.Beta];
             tspan = [0.0; 1.0];
             
             obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
index ca31dccc..3a64f2f6 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEParameters.m
@@ -1,7 +1,21 @@
-classdef AscherLinearDAEParameters
+classdef AscherLinearDAEParameters < otp.Parameters
     % Parameters for Ascher Linear DAE problem 
     properties
         % A scalar parameter $β$ in the linear model. It affects the stifness of the problem.
         Beta %MATLAB ONLY: (1,1) {mustBeNumeric} = 1
     end
+
+    methods
+        function obj = AscherLinearDAEParameters(varargin)
+            % Create a Ascher Linear DAE parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
+            %    value initializes that property.
+
+            obj = obj@otp.Parameters(varargin{:});
+        end
+    end
 end
diff --git a/toolbox/+otp/+cusp/+presets/Canonical.m b/toolbox/+otp/+cusp/+presets/Canonical.m
index 450cb77f..abe58344 100644
--- a/toolbox/+otp/+cusp/+presets/Canonical.m
+++ b/toolbox/+otp/+cusp/+presets/Canonical.m
@@ -8,7 +8,7 @@
     % b_i(0) &= 2 \sin\left( \frac{2 i \pi}{N} \right), \\
     % $$
     %
-    % for $i = 1, \dots, N$. The parameters are $\varepsilon = 10^{-4}$ and $\sigma = \frac{1}{144}$.
+    % for $i = 1, \dots, N$. The parameters are $ε = 10^{-4}$ and $σ = \frac{1}{144}$.
     
     methods
         function obj = Canonical(varargin)
@@ -20,8 +20,8 @@
             %    A variable number of name-value pairs. The accepted names are
             %
             %    - ``N`` – The number of cells in the spatial discretization.
-            %    - ``epsilon`` – Value of $\varepsilon$.
-            %    - ``sigma`` – Value of $\sigma$.
+            %    - ``epsilon`` – Value of $ε$.
+            %    - ``sigma`` – Value of $σ$.
             %
             % Returns
             % -------
@@ -30,22 +30,19 @@
 
             p = inputParser;
             p.addParameter('N', 32);
-            p.addParameter('epsilon', 1e-4);
-            p.addParameter('sigma', 1/144);
             p.parse(varargin{:});
-            opts = p.Results;
+            n = p.Results.N;
             
-            params = otp.cusp.CUSPParameters;
-            params.Epsilon = opts.epsilon;
-            params.Sigma = opts.sigma;
-            
-            ang = 2 * pi / opts.N * (1:opts.N).';
-            y0 = zeros(opts.N, 1);
-            a0 = -2*cos(ang);
-            b0 = 2*sin(ang);
+            ang = 2 * pi * (1:n).' / n;
+            y0 = zeros(n, 1);
+            a0 = -2 * cos(ang);
+            b0 = 2 * sin(ang);
 
             u0 = [y0; a0; b0];
             tspan = [0; 1.1];
+
+            unmatched = namedargs2cell(p.Unmatched);
+            params = otp.cusp.CUSPParameters('Epsilon', 1e-4, 'Sigma', 1/144, unmatched{:});
             
             obj = obj@otp.cusp.CUSPProblem(tspan, u0, params);
         end
diff --git a/toolbox/+otp/+cusp/CUSPParameters.m b/toolbox/+otp/+cusp/CUSPParameters.m
index 15b966e9..a56986c9 100644
--- a/toolbox/+otp/+cusp/CUSPParameters.m
+++ b/toolbox/+otp/+cusp/CUSPParameters.m
@@ -1,11 +1,25 @@
-classdef CUSPParameters
+classdef CUSPParameters < otp.Parameters
     % Parameters for the CUSP problem.
     properties
-        % The stiffness parameter $\varepsilon$ for the "cusp catastrophe" model.
+        % The stiffness parameter $ε$ for the "cusp catastrophe" model.
         Epsilon %MATLAB ONLY: (1,1) {mustBeNumeric}
 
-        % The diffusion coefficient $\sigma$ for all three variables.
+        % The diffusion coefficient $σ$ for all three variables.
         Sigma %MATLAB ONLY: (1,1) {mustBeNumeric}
     end
+
+    methods
+        function obj = CUSPParameters(varargin)
+            % Create a CUSP parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
+            %    value initializes that property.
+
+            obj = obj@otp.Parameters(varargin{:});
+        end
+    end
 end
 
diff --git a/toolbox/+otp/+cusp/CUSPProblem.m b/toolbox/+otp/+cusp/CUSPProblem.m
index 0794835c..dd86bbd0 100644
--- a/toolbox/+otp/+cusp/CUSPProblem.m
+++ b/toolbox/+otp/+cusp/CUSPProblem.m
@@ -4,10 +4,9 @@
     % The CUSP problem :cite:p:`HW96` (pp. 147-148) is the PDE
     %
     % $$
-    % \frac{\partial y}{\partial t} &= -\frac{1}{\varepsilon} (y^3 + a y + b)
-    % + \sigma \frac{\partial^2 y}{\partial x^2}, \\
-    % \frac{\partial a}{\partial t} &= b + 0.07 v + \sigma \frac{\partial^2 a}{\partial x^2}, \\
-    % \frac{\partial b}{\partial t} &= (1 - a^2) b - a - 0.4 y + 0.035 v + \sigma \frac{\partial^2 b}{\partial x^2},
+    % \frac{\partial y}{\partial t} &= -\frac{1}{ε} (y^3 + a y + b) + σ \frac{\partial^2 y}{\partial x^2}, \\
+    % \frac{\partial a}{\partial t} &= b + 0.07 v + σ \frac{\partial^2 a}{\partial x^2}, \\
+    % \frac{\partial b}{\partial t} &= (1 - a^2) b - a - 0.4 y + 0.035 v + σ \frac{\partial^2 b}{\partial x^2},
     % $$
     %
     %
@@ -15,22 +14,22 @@
     % conditions. Discretization with second order finite difference on a grid with $N$ cells gives
     %
     % $$
-    % y'_i &= -\frac{1}{\varepsilon} (y_i^3 + a_i y_i + b_i) + \sigma N^2 (y_{i-1} - 2 y_i + y_{i+1}) \\
-    % a'_i &= b_i + 0.07 v_i + \sigma N^2 (a_{i-1} - 2 a_i + a_{i+1}) \\
-    % b'_i &= (1 - a_i^2) b_i - a_i - 0.4 y_i + 0.035 v_i + \sigma N^2  (b_{i-1} - 2 b_i + b_{i+1}),
+    % y'_i &= -\frac{1}{ε} (y_i^3 + a_i y_i + b_i) + σ N^2 (y_{i-1} - 2 y_i + y_{i+1}) \\
+    % a'_i &= b_i + 0.07 v_i + σ N^2 (a_{i-1} - 2 a_i + a_{i+1}) \\
+    % b'_i &= (1 - a_i^2) b_i - a_i - 0.4 y_i + 0.035 v_i + σ N^2  (b_{i-1} - 2 b_i + b_{i+1}),
     % $$
     %
     % where $i = 1, \dots, N$. Values at cell indices $i=0, N+1$ are specificied by the periodic boundary conditions.
     %
     % Notes
     % -----
-    % +---------------------+----------------------------------------------------------------------------+
-    % | Type                | PDE                                                                        |
-    % +---------------------+----------------------------------------------------------------------------+
-    % | Number of Variables | arbitrary multiple of 3                                                    |
-    % +---------------------+----------------------------------------------------------------------------+
-    % | Stiff               | typically, depending on $\varepsilon$, $\sigma$, and number of grid points |
-    % +---------------------+----------------------------------------------------------------------------+
+    % +---------------------+-------------------------------------------------------------+
+    % | Type                | PDE                                                         |
+    % +---------------------+-------------------------------------------------------------+
+    % | Number of Variables | arbitrary multiple of 3                                     |
+    % +---------------------+-------------------------------------------------------------+
+    % | Stiff               | typically, depending on $ε$, $σ$, and number of grid points |
+    % +---------------------+-------------------------------------------------------------+
     %
     % Example
     % -------
@@ -44,7 +43,7 @@
     end
     
     properties (SetAccess = private)
-        % Right-hand side containing the diffusion terms and the reaction terms multiplied by $\varepsilon^{-1}$.
+        % Right-hand side containing the diffusion terms and the reaction terms multiplied by $ε^{-1}$.
         %
         % This partition of the RHS is used in :cite:p:`JM17`.
         %
@@ -53,7 +52,7 @@
         % RHSNonstiff
        RHSStiff
 
-        % Right-hand side containing the reaction terms not scaled by $\varepsilon^{-1}$.
+        % Right-hand side containing the reaction terms not scaled by $ε^{-1}$.
         %
         % This partition of the RHS is used in :cite:p:`JM17`.
         %
diff --git a/toolbox/+otp/+lorenz63/+presets/Canonical.m b/toolbox/+otp/+lorenz63/+presets/Canonical.m
index 4a8b36e4..3a87ac1f 100644
--- a/toolbox/+otp/+lorenz63/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz63/+presets/Canonical.m
@@ -15,23 +15,10 @@
             %    - ``sigma`` – Value of $σ$.
             %    - ``rho`` – Value of $ρ$.
             %    - ``beta`` – Value of $β$.
-            %
-
-            p = inputParser;
-            p.addParameter('sigma', 10);
-            p.addParameter('rho', 28);
-            p.addParameter('beta', 8/3);
-            p.parse(varargin{:});
-            opts = p.Results;
-
-            params = otp.lorenz63.Lorenz63Parameters;
-            
-            params.Sigma = opts.sigma;
-            params.Rho   = opts.rho;
-            params.Beta  = opts.beta;
             
             y0    = [0; 1; 0];
             tspan = [0 60];
+            params = otp.lorenz63.Lorenz63Parameters('Sigma', 10, 'Rho', 28, 'Beta', 8/3, varargin{:});
             
             obj = obj@otp.lorenz63.Lorenz63Problem(tspan, y0, params);            
         end
diff --git a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
index 3f525411..e7e75fd4 100644
--- a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
+++ b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
@@ -6,22 +6,13 @@
     methods
         function obj = LimitCycle
             % Create the LimitCycle Lorenz '63 problem object.
-            %
-
-            sigma = 10;
-            rho   = 350;
-            beta  = 8/3;
-    
-            params = otp.lorenz63.Lorenz63Parameters;
-            params.Sigma = sigma;
-            params.Rho   = rho;
-            params.Beta  = beta;
             
             % We use Lorenz's initial conditions and timespan as Strogatz
             % does not specify those in his book.
             
             y0    = [0; 1; 0];
             tspan = [0 60];
+            params = otp.lorenz63.Lorenz63Parameters('Sigma', 10, 'Rho', 350, 'Beta', 8/3);
             
             obj = obj@otp.lorenz63.Lorenz63Problem(tspan, y0, params);            
         end
diff --git a/toolbox/+otp/+lorenz63/+presets/Surprise.m b/toolbox/+otp/+lorenz63/+presets/Surprise.m
index 58f383ad..349b3e73 100644
--- a/toolbox/+otp/+lorenz63/+presets/Surprise.m
+++ b/toolbox/+otp/+lorenz63/+presets/Surprise.m
@@ -5,21 +5,11 @@
     methods
         function obj = Surprise
             % Create the Surprise Lorenz '63 problem object.
-            %
-
-            sigma = 10;
-            rho   = 100;
-            beta  = 8/3;
-
-            params = otp.lorenz63.Lorenz63Parameters;
-            params.Sigma = sigma;
-            params.Rho   = rho;
-            params.Beta  = beta;
             
             % Hand-picked initial conditions with the canonical timespan
-            
             y0    = [2; 1; 1];
             tspan = [0 60];
+            params = otp.lorenz63.Lorenz63Parameters('Sigma', 10, 'Rho', 100, 'Beta', 8/3);
             
             obj = obj@otp.lorenz63.Lorenz63Problem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+lorenz63/Lorenz63Parameters.m b/toolbox/+otp/+lorenz63/Lorenz63Parameters.m
index 919869ec..9f4a8aa6 100644
--- a/toolbox/+otp/+lorenz63/Lorenz63Parameters.m
+++ b/toolbox/+otp/+lorenz63/Lorenz63Parameters.m
@@ -1,4 +1,4 @@
-classdef Lorenz63Parameters
+classdef Lorenz63Parameters < otp.Parameters
     % Parameters for the Lorenz '63 problem.
 
     properties
@@ -11,4 +11,18 @@
         % A geometric factor.
         Beta %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
     end
+
+    methods
+        function obj = Lorenz63Parameters(varargin)
+            % Create a Lorenz '63 parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
+            %    value initializes that property.
+
+            obj = obj@otp.Parameters(varargin{:});
+        end
+    end
 end
diff --git a/toolbox/+otp/+lorenz96/+presets/Canonical.m b/toolbox/+otp/+lorenz96/+presets/Canonical.m
index a21c4784..c34522d5 100644
--- a/toolbox/+otp/+lorenz96/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz96/+presets/Canonical.m
@@ -12,31 +12,24 @@
             % varargin
             %    A variable number of name-value pairs. The accepted names are
             %
-            %    - ``Size`` – The size of the problem as a positive integer.
+            %    - ``N`` – The size of the problem as a positive integer.
             %    - ``Forcing`` – The forcing as a scalar, vector of N constants, or as a function.
             %
 
             p = inputParser;
-            p.addParameter('Size', 40, @isscalar);
-            p.addParameter('Forcing', 8);
-
+            p.addParameter('N', 40);
             p.parse(varargin{:});
-            
-            s = p.Results;
-            
-            n = s.Size;
-            
-            params = otp.lorenz96.Lorenz96Parameters;
-            params.F = s.Forcing;
+            n = p.Results.N;
             
             % We initialise the Lorenz96 model as in (Lorenz & Emanuel 1998)
-            
             y0 = 8*ones(n, 1);
-            
             y0(floor(n/2)) = 8.008;
             
             % roughly ten years in system time
             tspan = [0, 720];
+
+            unmatched = namedargs2cell(p.Unmatched);
+            params = otp.lorenz96.Lorenz96Parameters('F', 8, unmatched{:});
             
             obj = obj@otp.lorenz96.Lorenz96Problem(tspan, y0, params);
             
diff --git a/toolbox/+otp/+lorenz96/Lorenz96Parameters.m b/toolbox/+otp/+lorenz96/Lorenz96Parameters.m
index bb3ee4a6..4995b71c 100644
--- a/toolbox/+otp/+lorenz96/Lorenz96Parameters.m
+++ b/toolbox/+otp/+lorenz96/Lorenz96Parameters.m
@@ -1,8 +1,22 @@
-classdef Lorenz96Parameters
+classdef Lorenz96Parameters < otp.Parameters
     % Parameters for the Lorenz '96 problem.
     
     properties
         % A forcing function, scalar, or vector.
         F %MATLAB ONLY: {mustBeA(F, {'numeric','function_handle'})}
     end
+
+    methods
+        function obj = Lorenz96Parameters(varargin)
+            % Create a Lorenz '96 parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
+            %    value initializes that property.
+
+            obj = obj@otp.Parameters(varargin{:});
+        end
+    end
 end
diff --git a/toolbox/+otp/+oregonator/+presets/Canonical.m b/toolbox/+otp/+oregonator/+presets/Canonical.m
index 7c40b277..cad192bc 100644
--- a/toolbox/+otp/+oregonator/+presets/Canonical.m
+++ b/toolbox/+otp/+oregonator/+presets/Canonical.m
@@ -28,21 +28,9 @@
             % obj : Canonical
             %    The constructed problem.
 
-            p = inputParser;
-            p.addParameter('f', 1);
-            p.addParameter('q', 8.375e-6);
-            p.addParameter('s', 77.27);
-            p.addParameter('w', 0.1610);
-            p.parse(varargin{:});
-            opts = p.Results;
-
             tspan = [0, 360];
             y0 = [1; 2; 3];
-            params = otp.oregonator.OregonatorParameters;
-            params.F = 1;
-            params.Q = 8.375e-6;
-            params.S = 77.27;
-            params.W = 0.1610;
+            params = otp.oregonator.OregonatorParameters('F', 1, 'Q', 8.375e-6, 'S', 77.27, 'W', 0.1610, varargin{:});
             
             obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+oregonator/+presets/EnrightHull.m b/toolbox/+otp/+oregonator/+presets/EnrightHull.m
index 38eff49b..418ea2dd 100644
--- a/toolbox/+otp/+oregonator/+presets/EnrightHull.m
+++ b/toolbox/+otp/+oregonator/+presets/EnrightHull.m
@@ -23,11 +23,7 @@
 
             tspan = [0, 300];
             y0 = [4; 1.1; 4];
-            params = otp.oregonator.OregonatorParameters;
-            params.F = 1;
-            params.Q = 8.375e-6;
-            params.S = 77.27;
-            params.W = 0.1610;
+            params = otp.oregonator.OregonatorParameters('F', 1, 'Q', 8.375e-6, 'S', 77.27, 'W', 0.1610);
             
             obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m b/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
index 3dfc4c9b..eeae9466 100644
--- a/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
+++ b/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
@@ -25,11 +25,7 @@
 
             tspan = [0, 302.85805];
             y0 = [4; 1.331391; 2.852348];
-            params = otp.oregonator.OregonatorParameters;
-            params.F = 1;
-            params.Q = 8.375e-6;
-            params.S = 77.27;
-            params.W = 0.1610;
+            params = otp.oregonator.OregonatorParameters('F', 1, 'Q', 8.375e-6, 'S', 77.27, 'W', 0.1610);
             
             obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+oregonator/OregonatorParameters.m b/toolbox/+otp/+oregonator/OregonatorParameters.m
index ff354805..d5014b50 100644
--- a/toolbox/+otp/+oregonator/OregonatorParameters.m
+++ b/toolbox/+otp/+oregonator/OregonatorParameters.m
@@ -1,4 +1,4 @@
-classdef OregonatorParameters
+classdef OregonatorParameters < otp.Parameters
     % Parameters for the Oregonator problem.
     properties
         % The stoichiometric factor $f$.
@@ -13,5 +13,19 @@
         % The reaction constant $w$.
         W %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
     end
+
+    methods
+        function obj = OregonatorParameters(varargin)
+            % Create a Oregonator parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
+            %    value initializes that property.
+
+            obj = obj@otp.Parameters(varargin{:});
+        end
+    end
 end
 
diff --git a/toolbox/+otp/+protherorobinson/+presets/Canonical.m b/toolbox/+otp/+protherorobinson/+presets/Canonical.m
index ede364da..5e720645 100644
--- a/toolbox/+otp/+protherorobinson/+presets/Canonical.m
+++ b/toolbox/+otp/+protherorobinson/+presets/Canonical.m
@@ -1,44 +1,34 @@
 classdef Canonical < otp.protherorobinson.ProtheroRobinsonProblem
     % The Prothero–Robinson configuration from :cite:p:`PR74` (p. 159) based on :cite:p:`SLH70` (p. 272). It uses time
-    % span $t \in [0, 15]$, initial condition $y(0) = \phi(0)$, and parameters
+    % span $t \in [0, 15]$, initial condition $y(0) = φ(0)$, and parameters
     %
     % $$
-    % \lambda &= -200, \\
-    % \phi(t) &= 10 - (10 + t) e^{-t}.
+    % λ &= -200, \\
+    % φ(t) &= 10 - (10 + t) e^{-t}.
     % $$
 
     methods
-        function obj = Canonical(lambda, phi, dphi)
+        function obj = Canonical(varargin)
             % Create the Canonical Prothero–Robinson problem object.
             %
             % Parameters
             % ----------
-            % 
-            % lambda : numeric(1, 1), default=-200
-            %    The stiffness parameter and eigenvalue of the Jacobian $\lambda$.
-            % phi : function_handle, default=@(t)10-(10+t).*exp(-t)
-            %    The exact solution.
-            % dphi : function_handle, default=@(t)(9+t).*exp(-t)
-            %    The time derivative of $\phi(t)$.
+            % varargin
+            %    A variable number of name-value pairs. The accepted names are
+            %
+            %    - ``Lambda`` – The stiffness parameter and eigenvalue of the Jacobian $λ$.
+            %    - ``Phi`` – The exact solution.
+            %    - ``DPhi`` – The time derivative of $φ(t)$.
             %
             % Returns
             % -------
             % obj : ProtheroRobinsonProblem
             %    The constructed problem
-            if nargin < 1 || isempty(lambda)
-                lambda = -200;
-            end
-            if nargin < 2 || isempty(phi)
-                phi = @(t) 10 - (10 + t) .* exp(-t);
-                dphi = @(t) (9 + t) .* exp(-t);
-            end
             
-            params = otp.protherorobinson.ProtheroRobinsonParameters;
-            params.Lambda = lambda;
-            params.Phi = phi;
-            params.DPhi = dphi;
+            params = otp.protherorobinson.ProtheroRobinsonParameters('Lambda', -200, ...
+                'Phi', @(t) 10 - (10 + t) .* exp(-t), 'DPhi', @(t) (9 + t) .* exp(-t), varargin{:});
             tspan = [0, 15];
-            y0 = phi(tspan(1));
+            y0 = params.Phi(tspan(1));
             
             obj = obj@otp.protherorobinson.ProtheroRobinsonProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+protherorobinson/ProtheroRobinsonParameters.m b/toolbox/+otp/+protherorobinson/ProtheroRobinsonParameters.m
index 02c61931..25756a70 100644
--- a/toolbox/+otp/+protherorobinson/ProtheroRobinsonParameters.m
+++ b/toolbox/+otp/+protherorobinson/ProtheroRobinsonParameters.m
@@ -1,13 +1,27 @@
-classdef ProtheroRobinsonParameters
+classdef ProtheroRobinsonParameters < otp.Parameters
     % Parameters for the Prothero–Robinson problem.
     properties
-        % The stiffness parameter and eigenvalue of the Jacobian $\lambda$.
+        % The stiffness parameter and eigenvalue of the Jacobian $λ$.
         Lambda %MATLAB ONLY: (1,1) {mustBeFinite} = -1
 
-        % The function $\phi(t)$.
+        % The function $φ(t)$.
         Phi %MATLAB ONLY: {mustBeA(Phi, 'function_handle')} = @sin
 
-        % The time derivative of $\phi(t)$.
+        % The time derivative of $φ(t)$.
         DPhi %MATLAB ONLY: {mustBeA(DPhi, 'function_handle')} = @cos
     end
+
+    methods
+        function obj = ProtheroRobinsonParameters(varargin)
+            % Create a Prothero–Robinson parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
+            %    value initializes that property.
+
+            obj = obj@otp.Parameters(varargin{:});
+        end
+    end
 end
diff --git a/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m b/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m
index 22db1b5f..f3c65f7a 100644
--- a/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m
+++ b/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m
@@ -3,21 +3,21 @@
     %
     % The Prothero–Robinson problem :cite:p:`PR74` is the linear ODE
     % 
-    % $$y' = \lambda (y - \phi(t)) + \phi'(t).$$
+    % $$y' = λ (y - φ(t)) + φ'(t).$$
     %
     % This simple problem is used to test for order reduction and S-stability of time-stepping schemes. With initial
-    % condition $y(t_0) = \phi(t_0)$, the exact solution is $y(t) = \phi(t)$ independent of $\lambda$. Therefore, any
+    % condition $y(t_0) = φ(t_0)$, the exact solution is $y(t) = φ(t)$ independent of $λ$. Therefore, any
     % errors introduced by the numerical scheme can be measured easily.
     %
     % Notes
     % -----
-    % +---------------------+-----------------------------------+
-    % | Type                | ODE                               |
-    % +---------------------+-----------------------------------+
-    % | Number of Variables | 1                                 |
-    % +---------------------+-----------------------------------+
-    % | Stiff               | typically, depending on $\lambda$ |
-    % +---------------------+-----------------------------------+
+    % +---------------------+-----------------------------+
+    % | Type                | ODE                         |
+    % +---------------------+-----------------------------+
+    % | Number of Variables | 1                           |
+    % +---------------------+-----------------------------+
+    % | Stiff               | typically, depending on $λ$ |
+    % +---------------------+-----------------------------+
     %
     % Example
     % -------
diff --git a/toolbox/+otp/+robertson/+presets/Canonical.m b/toolbox/+otp/+robertson/+presets/Canonical.m
index bf921bc9..d769a73c 100644
--- a/toolbox/+otp/+robertson/+presets/Canonical.m
+++ b/toolbox/+otp/+robertson/+presets/Canonical.m
@@ -7,24 +7,26 @@
     % K_3 &= 10^4.
     % $$
     methods
-        function obj = Canonical
+        function obj = Canonical(varargin)
             % Create the Canonical Robertson problem object.
             %
             % Parameters
             % ----------
+            % varargin
+            %    A variable number of name-value pairs. The accepted names are
+            %
+            %    - ``K1`` – Value of $K_1$.
+            %    - ``K2`` – Value of $K_2$.
+            %    - ``K3`` – Value of $K_3$.
             %
             % Returns
             % -------
             % obj : RobertsonProblem
             %    The constructed problem.
             
-            params = otp.robertson.RobertsonParameters;
-            params.K1 = 0.04;
-            params.K2 = 3e7;
-            params.K3 = 1e4;
-            
             y0 = [1; 0; 0];
             tspan = [0, 40];
+            params = otp.robertson.RobertsonParameters('K1', 0.04, 'K2', 3e7, 'K3', 1e4, varargin{:});
             
             obj = obj@otp.robertson.RobertsonProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+robertson/+presets/ROBER.m b/toolbox/+otp/+robertson/+presets/ROBER.m
index 291da496..3df1ed30 100644
--- a/toolbox/+otp/+robertson/+presets/ROBER.m
+++ b/toolbox/+otp/+robertson/+presets/ROBER.m
@@ -6,21 +6,14 @@
         function obj = ROBER
             % Create the ROBER Robertson problem object.
             %
-            % Parameters
-            % ----------
-            %
             % Returns
             % -------
             % obj : RobertsonProblem
             %    The constructed problem.
             
-            params = otp.robertson.RobertsonParameters;
-            params.K1 = 0.04;
-            params.K2 = 3e7;
-            params.K3 = 1e4;
-            
             y0 = [1; 0; 0];
             tspan = [0; 1e11];
+            params = otp.robertson.RobertsonParameters('K1', 0.04, 'K2', 3e7, 'K3', 1e4);
             
             obj = obj@otp.robertson.RobertsonProblem(tspan, y0, params);
         end
diff --git a/toolbox/+otp/+robertson/RobertsonParameters.m b/toolbox/+otp/+robertson/RobertsonParameters.m
index a98905ef..fea738ea 100644
--- a/toolbox/+otp/+robertson/RobertsonParameters.m
+++ b/toolbox/+otp/+robertson/RobertsonParameters.m
@@ -1,4 +1,4 @@
-classdef RobertsonParameters
+classdef RobertsonParameters < otp.Parameters
     % Parameters for the Robertson problem.
     properties
         % The reaction rate $K_1$.
@@ -8,4 +8,18 @@
         % The reaction rate $K_3$.
         K3 %MATLAB ONLY: (1, 1) {mustBeReal, mustBeNonnegative}
     end
+
+    methods
+        function obj = RobertsonParameters(varargin)
+            % Create a Robertson parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
+            %    value initializes that property.
+
+            obj = obj@otp.Parameters(varargin{:});
+        end
+    end
 end
diff --git a/toolbox/+otp/Parameters.m b/toolbox/+otp/Parameters.m
new file mode 100644
index 00000000..6669afc2
--- /dev/null
+++ b/toolbox/+otp/Parameters.m
@@ -0,0 +1,28 @@
+classdef Parameters
+    % A superclass for the parameters of a differential equation.
+    %
+    % Subclasses of this class specify mutable properties associated with an :class:`otp.Problem`.
+    %
+    % See Also
+    % --------
+    % otp.Problem.Parameters
+
+    methods
+        function obj = Parameters(varargin)
+            % Create a parameters object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. A name can be any property defined by a subclass, and the
+            %    subsequent value initializes that property.
+            extras = struct(varargin{:});
+            fields = fieldnames(extras);
+            
+            for i = 1:length(fields)
+                f = fields{i};
+                obj.(f) = extras.(f);
+            end
+        end
+    end
+end
diff --git a/toolbox/+otp/Problem.m b/toolbox/+otp/Problem.m
index 6168facc..19e8d6bb 100644
--- a/toolbox/+otp/Problem.m
+++ b/toolbox/+otp/Problem.m
@@ -34,7 +34,7 @@
         Y0 %MATLAB ONLY: (:,1) {otp.utils.validation.mustBeNumerical}
         
         % Additional variables to pass to the F function
-        Parameters %MATLAB ONLY: (1,1)
+        Parameters %MATLAB ONLY: (1,1) otp.Parameters
         
         % The dimension of the ODE
         NumVars
diff --git a/toolbox/+otp/RHS.m b/toolbox/+otp/RHS.m
index 15e4d492..869bf940 100644
--- a/toolbox/+otp/RHS.m
+++ b/toolbox/+otp/RHS.m
@@ -291,7 +291,7 @@
             % F : function_handle
             %    The function handle for $f$ in the differential equation $M(t, y) y' = f(t, y)$.
             % varargin
-            %    A variable number of name-value pairs. A name can be any property of :class:`RHS`, and the subsequent
+            %    A variable number of name-value pairs. A name can be any property of this class, and the subsequent
             %    value initializes that property.
             %
             % Warning

From e6ead75f9e4b86c6c8d3be8577044fbe29a7fb00 Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Sun, 5 Nov 2023 21:59:51 -0800
Subject: [PATCH 16/20] Doc consistency fixes

---
 docs/conf.py                                  |  1 +
 docs/references.bib                           |  9 ++--
 .../+ascherlineardae/+presets/Canonical.m     |  7 ++-
 .../+otp/+ascherlineardae/+presets/Petzold.m  |  9 ++--
 .../+otp/+ascherlineardae/+presets/Stiff.m    |  9 ++--
 .../+ascherlineardae/AscherLinearDAEProblem.m | 11 +++--
 toolbox/+otp/+cusp/+presets/Canonical.m       | 15 ++----
 toolbox/+otp/+cusp/CUSPProblem.m              | 10 ++--
 toolbox/+otp/+lorenz63/+presets/Canonical.m   | 10 ++--
 toolbox/+otp/+lorenz63/+presets/LimitCycle.m  |  7 +--
 toolbox/+otp/+lorenz63/+presets/Surprise.m    |  5 +-
 toolbox/+otp/+lorenz63/Lorenz63Problem.m      | 49 +++++++++----------
 toolbox/+otp/+lorenz96/+presets/Canonical.m   |  9 ++--
 toolbox/+otp/+lorenz96/+presets/PopovSandu.m  | 41 ++++++----------
 toolbox/+otp/+lorenz96/Lorenz96Problem.m      | 30 ++++++------
 toolbox/+otp/+oregonator/+presets/Canonical.m | 12 ++---
 .../+otp/+oregonator/+presets/EnrightHull.m   | 14 ++----
 .../+oregonator/+presets/GottwaldWanner.m     | 14 ++----
 toolbox/+otp/+oregonator/OregonatorProblem.m  |  6 +--
 .../+protherorobinson/+presets/Canonical.m    |  9 +---
 .../ProtheroRobinsonProblem.m                 |  6 +--
 toolbox/+otp/+robertson/+presets/Canonical.m  | 22 +++------
 toolbox/+otp/+robertson/+presets/ROBER.m      | 10 +---
 toolbox/+otp/+robertson/RobertsonProblem.m    | 40 +++++++--------
 toolbox/+otp/RHS.m                            |  2 +
 25 files changed, 133 insertions(+), 224 deletions(-)

diff --git a/docs/conf.py b/docs/conf.py
index c5f19942..2ab3a739 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -26,6 +26,7 @@
 primary_domain = 'mat'
 matlab_src_dir = '../toolbox'
 matlab_keep_package_prefix = False
+matlab_auto_link = 'basic'
 
 autodoc_default_options = {
     'member-order': 'bysource',
diff --git a/docs/references.bib b/docs/references.bib
index 2a00c9cb..728b7b8a 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -55,17 +55,16 @@ @Article{PR74
 }
 
 
-@Article{Rob67,
+@Article{Rob66,
     title={The Solution of a Set of Reaction Rate Equations},
     author={Robertson, H.H.},
     journal={Walsh, J., Ed., Numerical Analysis, An introduction},
-    year={1967},
+    year={1966},
     pages= {178--182},
-    publisher = "Academic Press",
-    address   = "Cambridge, Massachusetts",
+    volume={178182},
+    publisher = {Academic Press}
 }
 
-
 @article{SLH70,
     author = {Seinfeld, J. H. and Lapidus, Leon and Hwang, Myungkyu},
     title = {Review of Numerical Integration Techniques for Stiff Ordinary Differential Equations},
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
index 1f6a35c5..d6f697b4 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Canonical.m
@@ -1,10 +1,10 @@
 classdef Canonical < otp.ascherlineardae.AscherLinearDAEProblem
-    % The problem defined by Uri Ascher in :cite:p:`Asc89` (sec. 2) which uses time span $t \in [0, 1]$ and intial 
+    % The problem defined by Uri Ascher in :cite:p:`Asc89` (sec. 2) which uses time span $t ∈ [0, 1]$ and intial 
     % condition $[y_0, z_0]^T = [1, β]^T$.
-    % 
+
     methods
         function obj = Canonical(varargin)
-            % Create the Canonical CUSP problem object.
+            % Create the Canonical Ascher Linear DAE problem object.
             %
             % Parameters
             % ----------
@@ -16,7 +16,6 @@
             params = otp.ascherlineardae.AscherLinearDAEParameters('Beta', 1, varargin{:});
             y0 = [1; params.Beta];
             tspan = [0.0; 1.0];
-            
             obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
         end
     end
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
index 0e67717e..14e9f89b 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Petzold.m
@@ -1,15 +1,14 @@
 classdef Petzold < otp.ascherlineardae.AscherLinearDAEProblem
     % The Petzold DAE example :cite:p:`Pet86` as a special case of the Ascher linear DAE problem. This preset uses time
-    % span $t \in [0, 1]$ and $β = 0 $ with the initial condition $[y_0, z_0]^T = [1, 0]^T $.
-    % 
+    % span $t ∈ [0, 1]$ and $β = 0 $ with the initial condition $[y_0, z_0]^T = [1, 0]^T$.
+
     methods
-        function obj = Petzold
-            % Create the Petzold example of the Ascher linear DAE problem object.
+        function obj = Petzold()
+            % Create the Petzold Ascher Linear DAE problem object.
             
             params = otp.ascherlineardae.AscherLinearDAEParameters('Beta', 0);
             y0 = [1; params.Beta];
             tspan = [0.0; 1.0];
-            
             obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
         end
     end
diff --git a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
index a32cbfe7..5e102971 100644
--- a/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
+++ b/toolbox/+otp/+ascherlineardae/+presets/Stiff.m
@@ -1,15 +1,14 @@
 classdef Stiff < otp.ascherlineardae.AscherLinearDAEProblem
     % The Stiff example from :cite:p:`Asc89`. A variant of the Ascher linear DAE problem which uses time span 
-    % $t \in [0, 1]$ and $β = 100$ with the initial condition $[y_0, z_0]^T = [1, 100]^T$.
-    % 
+    % $t ∈ [0, 1]$ and $β = 100$ with the initial condition $[y_0, z_0]^T = [1, 100]^T$.
+
     methods
-        function obj = Stiff
-            % Create the stiff example of the Ascher linear DAE problem object.
+        function obj = Stiff()
+            % Create the stiff Ascher Linear DAE problem object.
 
             params = otp.ascherlineardae.AscherLinearDAEParameters('Beta', 100);
             y0 = [1; params.Beta];
             tspan = [0.0; 1.0];
-            
             obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
         end
     end
diff --git a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
index e68c84af..7ceea46b 100644
--- a/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
+++ b/toolbox/+otp/+ascherlineardae/AscherLinearDAEProblem.m
@@ -1,5 +1,5 @@
 classdef AscherLinearDAEProblem < otp.Problem
-    % A linear differential-algebraic problem with time-dependant mass matrix.
+    % A linear differential-algebraic problem with a time-dependant mass matrix.
     %
     % The Ascher linear DAE Problem :cite:p:`Asc89` is an index-1 differential-agebraic equation given by
     %
@@ -7,10 +7,10 @@
     % \begin{bmatrix}
     % 1 & -t \\
     % 0 & 0
-    % \end{bmatrix} \begin{bmatrix} y'(t) \\ z'(t) \end{bmatrix} = \left[ \begin{array}{cc}
+    % \end{bmatrix} \begin{bmatrix} y'(t) \\ z'(t) \end{bmatrix} = \begin{bmatrix}
     % -1 & 1+t \\
     % β & -1-β t
-    % \end{array}\right] \begin{bmatrix} y(t) \\ z(t) \end{bmatrix} + \begin{bmatrix}
+    % \end{bmatrix} \begin{bmatrix} y(t) \\ z(t) \end{bmatrix} + \begin{bmatrix}
     % 0 \\
     % \sin(t)
     % \end{bmatrix}.
@@ -20,10 +20,11 @@
     %
     % $$
     % \begin{bmatrix} y(t)\\ z(t) \end{bmatrix} = \begin{bmatrix}
-    % t \sin(t) + (1 + β  t) e^{-t}\\
-    % β  e^{-t} + \sin(t)
+    % t \sin(t) + (1 + β t) e^{-t}\\
+    % β e^{-t} + \sin(t)
     % \end{bmatrix}.
     % $$
+    %
     % This DAE problem can be used to investigate the convergence of implcit time-stepping methods due to its stiffness
     % and time-dependant mass matrix.
     %
diff --git a/toolbox/+otp/+cusp/+presets/Canonical.m b/toolbox/+otp/+cusp/+presets/Canonical.m
index abe58344..01900d0b 100644
--- a/toolbox/+otp/+cusp/+presets/Canonical.m
+++ b/toolbox/+otp/+cusp/+presets/Canonical.m
@@ -1,14 +1,14 @@
 classdef Canonical < otp.cusp.CUSPProblem
-    % The CUSP configuration from :cite:p:`HW96` (pp. 147-148) which uses time span $t \in [0, 1.1]$, $N = 32$ grid
-    % cells, and initial conditions
+    % The CUSP configuration from :cite:p:`HW96` (pp. 147-148) which uses time span $t ∈ [0, 1.1]$, $N = 32$ grid cells,
+    % and initial conditions
     %
     % $$
     % y_i(0) &= 0 ,\\
-    % a_i(0) &= -2 \cos\left( \frac{2 i \pi}{N} \right), \\
-    % b_i(0) &= 2 \sin\left( \frac{2 i \pi}{N} \right), \\
+    % a_i(0) &= -2 \cos\left( \frac{2 i π}{N} \right), \\
+    % b_i(0) &= 2 \sin\left( \frac{2 i π}{N} \right), \\
     % $$
     %
-    % for $i = 1, \dots, N$. The parameters are $ε = 10^{-4}$ and $σ = \frac{1}{144}$.
+    % for $i = 1, …, N$. The parameters are $ε = 10^{-4}$ and $σ = \frac{1}{144}$.
     
     methods
         function obj = Canonical(varargin)
@@ -22,11 +22,6 @@
             %    - ``N`` – The number of cells in the spatial discretization.
             %    - ``epsilon`` – Value of $ε$.
             %    - ``sigma`` – Value of $σ$.
-            %
-            % Returns
-            % -------
-            % obj : CUSPProblem
-            %    The constructed problem.
 
             p = inputParser;
             p.addParameter('N', 32);
diff --git a/toolbox/+otp/+cusp/CUSPProblem.m b/toolbox/+otp/+cusp/CUSPProblem.m
index dd86bbd0..d637351f 100644
--- a/toolbox/+otp/+cusp/CUSPProblem.m
+++ b/toolbox/+otp/+cusp/CUSPProblem.m
@@ -10,7 +10,7 @@
     % $$
     %
     %
-    % where $v = u / (u + 0.1)$ and $u = (y - 0.7)(y - 1.3)$. The spatial domain $x \in [0, 1]$ has period boundary
+    % where $v = u / (u + 0.1)$ and $u = (y - 0.7)(y - 1.3)$. The spatial domain $x ∈ [0, 1]$ has period boundary
     % conditions. Discretization with second order finite difference on a grid with $N$ cells gives
     %
     % $$
@@ -19,7 +19,7 @@
     % b'_i &= (1 - a_i^2) b_i - a_i - 0.4 y_i + 0.035 v_i + σ N^2  (b_{i-1} - 2 b_i + b_{i+1}),
     % $$
     %
-    % where $i = 1, \dots, N$. Values at cell indices $i=0, N+1$ are specificied by the periodic boundary conditions.
+    % where $i = 1, …, N$. Values at cell indices $i=0, N+1$ are specificied by the periodic boundary conditions.
     %
     % Notes
     % -----
@@ -88,11 +88,7 @@
             %    The initial conditions.
             % parameters : CUSPParameters
             %    The parameters.
-            %
-            % Returns
-            % -------
-            % obj : CUSPProblem
-            %    The constructed problem.
+            
             obj@otp.Problem('CUSP', [], timeSpan, y0, parameters);
         end
     end
diff --git a/toolbox/+otp/+lorenz63/+presets/Canonical.m b/toolbox/+otp/+lorenz63/+presets/Canonical.m
index 3a87ac1f..ca63c3e8 100644
--- a/toolbox/+otp/+lorenz63/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz63/+presets/Canonical.m
@@ -1,6 +1,5 @@
 classdef Canonical < otp.lorenz63.Lorenz63Problem
-    % Original Lorenz '63 preset presented in :cite:p:`Lor63`
-    % which uses time span $t \in [0, 60]$, $σ = 10$, $ρ = 28$, 
+    % Original Lorenz '63 preset presented in :cite:p:`Lor63` which uses time span $t ∈ [0, 60]$, $σ = 10$, $ρ = 28$, 
     % $β = 8/3$, and intial conditions $y_0 = [0, 1, 0]^T$.
 
     methods
@@ -12,14 +11,13 @@
             % varargin
             %    A variable number of name-value pairs. The accepted names are
             %
-            %    - ``sigma`` – Value of $σ$.
-            %    - ``rho`` – Value of $ρ$.
-            %    - ``beta`` – Value of $β$.
+            %    - ``Sigma`` – Value of $σ$.
+            %    - ``Rho`` – Value of $ρ$.
+            %    - ``Beta`` – Value of $β$.
             
             y0    = [0; 1; 0];
             tspan = [0 60];
             params = otp.lorenz63.Lorenz63Parameters('Sigma', 10, 'Rho', 28, 'Beta', 8/3, varargin{:});
-            
             obj = obj@otp.lorenz63.Lorenz63Problem(tspan, y0, params);            
         end
     end
diff --git a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
index e7e75fd4..9981ba47 100644
--- a/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
+++ b/toolbox/+otp/+lorenz63/+presets/LimitCycle.m
@@ -1,7 +1,6 @@
 classdef LimitCycle < otp.lorenz63.Lorenz63Problem
-    % Lorenz '63 preset limit cycle, a non-chaotic preset, from :cite:p:`Str18` 
-    % which uses time span $t \in [0, 60]$, $σ = 10$, $ρ = 350$, 
-    % $β = 8/3$, and intial conditions $y_0 = [0, 1, 0]^T$.
+    % Lorenz '63 preset limit cycle, a non-chaotic preset, from :cite:p:`Str18` which uses time span $t ∈ [0, 60]$,
+    % $σ = 10$, $ρ = 350$, $β = 8/3$, and intial conditions $y_0 = [0, 1, 0]^T$.
     
     methods
         function obj = LimitCycle
@@ -9,11 +8,9 @@
             
             % We use Lorenz's initial conditions and timespan as Strogatz
             % does not specify those in his book.
-            
             y0    = [0; 1; 0];
             tspan = [0 60];
             params = otp.lorenz63.Lorenz63Parameters('Sigma', 10, 'Rho', 350, 'Beta', 8/3);
-            
             obj = obj@otp.lorenz63.Lorenz63Problem(tspan, y0, params);            
         end
     end
diff --git a/toolbox/+otp/+lorenz63/+presets/Surprise.m b/toolbox/+otp/+lorenz63/+presets/Surprise.m
index 349b3e73..40236a0c 100644
--- a/toolbox/+otp/+lorenz63/+presets/Surprise.m
+++ b/toolbox/+otp/+lorenz63/+presets/Surprise.m
@@ -1,6 +1,6 @@
 classdef Surprise < otp.lorenz63.Lorenz63Problem
-    % Lorenz '63 preset 'surprise' from :cite:p:`Str18` which uses time span $t \in [0, 60]$, 
-    % $σ = 10$, $ρ = 100$, $β = 8/3$, and intial conditions $y_0 = [2, 1, 1]^T$.
+    % Lorenz '63 preset 'surprise' from :cite:p:`Str18` which uses time span $t ∈ [0, 60]$,  $σ = 10$, $ρ = 100$,
+    % $β = 8/3$, and intial conditions $y_0 = [2, 1, 1]^T$.
 
     methods
         function obj = Surprise
@@ -10,7 +10,6 @@
             y0    = [2; 1; 1];
             tspan = [0 60];
             params = otp.lorenz63.Lorenz63Parameters('Sigma', 10, 'Rho', 100, 'Beta', 8/3);
-            
             obj = obj@otp.lorenz63.Lorenz63Problem(tspan, y0, params);
         end
     end
diff --git a/toolbox/+otp/+lorenz63/Lorenz63Problem.m b/toolbox/+otp/+lorenz63/Lorenz63Problem.m
index a9c41122..ebcf8b28 100644
--- a/toolbox/+otp/+lorenz63/Lorenz63Problem.m
+++ b/toolbox/+otp/+lorenz63/Lorenz63Problem.m
@@ -9,22 +9,19 @@
     % z' &= xy - βz,
     % $$
     %
-    % exhibits chaotic behavior for certain values of the parameters.
-    % Here $x$ roughly corresponds to the rate of convection of a fluid,
-    % $y$ corresponds to temperature variation in one direction,
-    % and $z$ is temerature variation in the other direction.
-    %
-    % For a full problem description see :cite:p:`Str18`.
+    % exhibits chaotic behavior for certain values of the parameters. Here $x$ roughly corresponds to the rate of
+    % convection of a fluid, $y$ corresponds to temperature variation in one direction, and $z$ is temerature variation
+    % in the other direction. For a full problem description see :cite:p:`Str18`.
     %
     % Notes
     % -----
-    % +---------------------+-----------------------------------------------------------+
-    % | Type                | ODE                                                       |
-    % +---------------------+-----------------------------------------------------------+
-    % | Number of Variables | 3                                                         |
-    % +---------------------+-----------------------------------------------------------+
-    % | Stiff               | not typically, depending on $σ$, $ρ$, and $β$             |
-    % +---------------------+-----------------------------------------------------------+
+    % +---------------------+-----------------------------------------------+
+    % | Type                | ODE                                           |
+    % +---------------------+-----------------------------------------------+
+    % | Number of Variables | 3                                             |
+    % +---------------------+-----------------------------------------------+
+    % | Stiff               | not typically, depending on $σ$, $ρ$, and $β$ |
+    % +---------------------+-----------------------------------------------+
     %
     % Example
     % -------
@@ -32,9 +29,9 @@
     % >>> sol = problem.solve('MaxStep', 1e-3);
     % >>> problem.plotPhaseSpace(sol);
     %
-    % See also
+    % See Also
     % --------
-    % :doc:`Lorenz '96 Problem </problems/Lorenz96Problem>`
+    % otp.lorenz96.Lorenz96Problem
     
     methods
         function obj = Lorenz63Problem(timeSpan, y0, parameters)
@@ -48,11 +45,7 @@
             %    The initial conditions.
             % parameters : Lorenz63Parameters
             %    The parameters.
-            %
-            % Returns
-            % -------
-            % obj : Lorenz63Problem
-            %    The constructed problem.
+
             obj@otp.Problem('Lorenz Equations', 3, timeSpan, y0, parameters);
         end
     end
@@ -64,12 +57,16 @@ function onSettingsChanged(obj)
             beta  = obj.Parameters.Beta;
             
             obj.RHS = otp.RHS(@(t, y) otp.lorenz63.f(t, y, sigma, rho, beta), ...
-                'Jacobian',                     @(t, y) otp.lorenz63.jacobian(t, y, sigma, rho, beta), ...
-                'JacobianVectorProduct',        @(t, y, v) otp.lorenz63.jacobianVectorProduct(t, y, v, sigma, rho, beta), ...
-                'JacobianAdjointVectorProduct', @(t, y, v) otp.lorenz63.jacobianAdjointVectorProduct(t, y, v, sigma, rho, beta), ...
-                'PartialDerivativeParameters',  @(t, y) otp.lorenz63.partialDerivativeParameters(t, y, sigma, rho, beta), ...
-                'HessianVectorProduct',         @(t, y, u, v) otp.lorenz63.hessianVectorProduct(t, y, u, v, sigma, rho, beta), ...
-                'HessianAdjointVectorProduct',  @(t, y, u, v) otp.lorenz63.hessianAdjointVectorProduct(t, y, u, v, sigma, rho, beta), ...
+                'Jacobian', @(t, y) otp.lorenz63.jacobian(t, y, sigma, rho, beta), ...
+                'JacobianVectorProduct', @(t, y, v) otp.lorenz63.jacobianVectorProduct(t, y, v, sigma, rho, beta), ...
+                'JacobianAdjointVectorProduct', ...
+                @(t, y, v) otp.lorenz63.jacobianAdjointVectorProduct(t, y, v, sigma, rho, beta), ...
+                'PartialDerivativeParameters', ...
+                @(t, y) otp.lorenz63.partialDerivativeParameters(t, y, sigma, rho, beta), ...
+                'HessianVectorProduct', ...
+                @(t, y, u, v) otp.lorenz63.hessianVectorProduct(t, y, u, v, sigma, rho, beta), ...
+                'HessianAdjointVectorProduct', ...
+                @(t, y, u, v) otp.lorenz63.hessianAdjointVectorProduct(t, y, u, v, sigma, rho, beta), ...
                 'Vectorized', 'on');
         end
         
diff --git a/toolbox/+otp/+lorenz96/+presets/Canonical.m b/toolbox/+otp/+lorenz96/+presets/Canonical.m
index c34522d5..d1094e36 100644
--- a/toolbox/+otp/+lorenz96/+presets/Canonical.m
+++ b/toolbox/+otp/+lorenz96/+presets/Canonical.m
@@ -1,11 +1,10 @@
 classdef Canonical < otp.lorenz96.Lorenz96Problem
-    % Original Lorenz '96 preset presented in :cite:p:`Lor96`
-    % which uses time span $t \in [0, 720]$, $N = 40$, $F=8$, and initial
-    % conditions of $y_i = 8$ for all $i$ except for $y_{\lfloor N/2 \rfloor}=8.008$. 
+    % Original Lorenz '96 preset presented in :cite:p:`Lor96` which uses time span $t ∈ [0, 720]$, $N = 40$, $F=8$, and
+    % initial conditions of $y_i = 8$ for all $i$ except for $y_{⌊N/2⌋}=8.008$. 
 
     methods
         function obj = Canonical(varargin)
-            % Create a Canonial problem object.
+            % Create the Canonical Lorenz '96 problem object.
             %
             % Parameters
             % ----------
@@ -14,7 +13,6 @@
             %
             %    - ``N`` – The size of the problem as a positive integer.
             %    - ``Forcing`` – The forcing as a scalar, vector of N constants, or as a function.
-            %
 
             p = inputParser;
             p.addParameter('N', 40);
@@ -32,7 +30,6 @@
             params = otp.lorenz96.Lorenz96Parameters('F', 8, unmatched{:});
             
             obj = obj@otp.lorenz96.Lorenz96Problem(tspan, y0, params);
-            
         end
     end
 end
diff --git a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
index 2b7b4f4f..b1beb936 100644
--- a/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
+++ b/toolbox/+otp/+lorenz96/+presets/PopovSandu.m
@@ -1,60 +1,47 @@
 classdef PopovSandu < otp.lorenz96.Lorenz96Problem
-    % A preset that has a cyclic forcing function that is different for
-    % every variable. This preset was created for :cite:p:`PS19`.
-    % This preset uses time span $t \in [0, 720]$, $N = 40$, and initial
-    % conditions of $y_i = 8$ for all $i$ except for $y_{\lfloor N/2 \rfloor}=8.008$. 
-    % The forcing, as a function of time is given by
+    % A preset that has a cyclic forcing function that is different for every variable. This preset was created for
+    % :cite:p:`PS19`. This preset uses time span $t ∈ [0, 720]$, $N = 40$, and initial conditions of $y_i = 8$ for all
+    % $i$ except for $y_{⌊N/2⌋}=8.008$. The forcing as a function of time is given by
     %
     % $$
-    % f(t) = 8 + 4\cos(2 \pi \omega (t+\text{mod}(i - 1, q)/q))
+    % f(t) = 8 + 4\cos(2 π ω (t+\text{mod}(i - 1, q)/q))
     % $$
     %
-    % where by default $q=4$ is the number of partitions, and $\omega = 1$
-    % is a forcing period of one time unit.
+    % where by default $q=4$ is the number of partitions, and $ω = 1$ is a forcing period of one time unit.
     
     methods
         function obj = PopovSandu(varargin)
-            % Create a PopovSandu problem object.
+            % Create the PopovSandu Lorenz '96 problem object.
             %
             % Parameters
             % ----------
             % varargin
             %    A variable number of name-value pairs. The accepted names are
             %
-            %    - ``Size`` – The size of the problem as a positive integer.
+            %    - ``N`` – The size of the problem as a positive integer.
             %    - ``Partitions`` – The number of partitions into which to divide the variables.
             %    - ``ForcingPeriod`` – The period of the forcing function in radians per unit time.
-            %
             
             p = inputParser;
-            
-            p.addParameter('Size', 40, @isscalar);
-            p.addParameter('Partitions', 4, @isscalar);
-            p.addParameter('ForcingPeriod', 1, @isscalar); % default corresponds to 5 days
-
+            p.addParameter('N', 40);
+            p.addParameter('Partitions', 4);
+            p.addParameter('ForcingPeriod', 1); % default corresponds to 5 days
             p.parse(varargin{:});
             
-            s = p.Results;
-            
-            n = s.Size;
-            q = s.Partitions;
-            omega = s.ForcingPeriod;
-           
-            F = @(t) 8 + 4*cos(2*pi*omega*(t + mod((1:n) - 1, q)/q)).';
+            n = p.Results.N;
+            q = p.Results.Partitions;
+            omega = p.Results.ForcingPeriod;
             
-            params = otp.lorenz96.Lorenz96Parameters;
-            params.F = F;            
+            params = otp.lorenz96.Lorenz96Parameters('F', @(t) 8 + 4*cos(2*pi*omega*(t + mod((1:n) - 1, q)/q)).');
             
             % We initialise the Lorenz96 model as in (Lorenz & Emanuel 1998)
             
             y0 = 8*ones(n, 1);
-            
             y0(floor(n/2)) = 8.008;
             
             tspan = [0, 720]; % 3600 days
             
             obj = obj@otp.lorenz96.Lorenz96Problem(tspan, y0, params);
-            
         end
     end
 end
diff --git a/toolbox/+otp/+lorenz96/Lorenz96Problem.m b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
index 459e9de8..694c96db 100644
--- a/toolbox/+otp/+lorenz96/Lorenz96Problem.m
+++ b/toolbox/+otp/+lorenz96/Lorenz96Problem.m
@@ -1,26 +1,24 @@
 classdef Lorenz96Problem <  otp.Problem
-    % A chaotic system modeling nonlinear transfer of a dimensionless
-    % quantity along a cyclic one dimensional domain.
+    % A chaotic system modeling nonlinear transfer of a dimensionless quantity along a cyclic one dimensional domain.
     % 
     % The $N$ variable dynamics :cite:p:`Lor96` are represented by the equation,
     %
     % $$
-    % y_i' = -y_{i-1}\left(y_{i-2} - y_{i+1}\right) - y_i + f(t), \quad i = 1, \dots, N,
+    % y_i' = -y_{i-1} (y_{i-2} - y_{i+1}) - y_i + f(t), \qquad i = 1, …, N,
     % $$
     % 
-    % where $y_0 = y_N$, $y_{-1} = y_{N - 1}$, and $y_{N + 1} = y_2$, exhibits 
-    % chaotic behavior for certain pairs of values of the dimension $N$ and
-    % forcing function $f$.
+    % where $y_0 = y_N$, $y_{-1} = y_{N - 1}$, and $y_{N + 1} = y_2$, exhibits chaotic behavior for certain pairs of
+    % values of the dimension $N$ and forcing function $f$.
     %
     % Notes
     % -----
-    % +---------------------+-----------------------------------------------------------+
-    % | Type                | ODE                                                       |
-    % +---------------------+-----------------------------------------------------------+
-    % | Number of Variables | $N$ for any postive integer four or greater               |
-    % +---------------------+-----------------------------------------------------------+
-    % | Stiff               | no                                                        |
-    % +---------------------+-----------------------------------------------------------+
+    % +---------------------+---------------------------------------------+
+    % | Type                | ODE                                         |
+    % +---------------------+---------------------------------------------+
+    % | Number of Variables | $N$ for any postive integer four or greater |
+    % +---------------------+---------------------------------------------+
+    % | Stiff               | no                                          |
+    % +---------------------+---------------------------------------------+
     %
     % Example
     % -------
@@ -28,9 +26,9 @@
     % >>> sol = problem.solve();
     % >>> problem.movie(sol);
     %
-    % See also
+    % See Also
     % --------
-    % :doc:`Lorenz '63 Problem </problems/Lorenz63Problem>`
+    % otp.lorenz63.Lorenz63Problem
     
     properties (SetAccess = private)
         DistanceFunction
@@ -48,7 +46,7 @@
             %    The initial conditions.
             % parameters : Lorenz96Parameters
             %    The parameters.
-            %
+
             obj@otp.Problem('Lorenz 96', [], timeSpan, y0, parameters);
         end
     end
diff --git a/toolbox/+otp/+oregonator/+presets/Canonical.m b/toolbox/+otp/+oregonator/+presets/Canonical.m
index cad192bc..fbbb273e 100644
--- a/toolbox/+otp/+oregonator/+presets/Canonical.m
+++ b/toolbox/+otp/+oregonator/+presets/Canonical.m
@@ -1,17 +1,17 @@
 classdef Canonical < otp.oregonator.OregonatorProblem
-    % The Oregonator configuration used in :cite:p:`HW96` (p. 144) called OREGO. It uses time span $t \in [0, 360]$,
+    % The Oregonator configuration used in :cite:p:`HW96` (p. 144) called OREGO. It uses time span $t ∈ [0, 360]$,
     % initial condition $y_0 = [1, 2, 3]^T$, and parameters
     %
     % $$
     % f &= 1, \\
-    % q &= 8.375 \cdot 10^{-6}, \\
+    % q &= 8.375 \times 10^{-6}, \\
     % s &= 77.27, \\
     % w &= 0.1610.
     % $$
 
     methods
         function obj = Canonical(varargin)
-            % Create the Canonical problem object.
+            % Create the Canonical Oregonator problem object.
             %
             % Parameters
             % ----------
@@ -22,16 +22,10 @@
             %    - ``q`` – Value of $q$.
             %    - ``s`` – Value of $s$.
             %    - ``w`` – Value of $w$.
-            %
-            % Returns
-            % -------
-            % obj : Canonical
-            %    The constructed problem.
 
             tspan = [0, 360];
             y0 = [1; 2; 3];
             params = otp.oregonator.OregonatorParameters('F', 1, 'Q', 8.375e-6, 'S', 77.27, 'W', 0.1610, varargin{:});
-            
             obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
         end
     end
diff --git a/toolbox/+otp/+oregonator/+presets/EnrightHull.m b/toolbox/+otp/+oregonator/+presets/EnrightHull.m
index 418ea2dd..c8a4cdde 100644
--- a/toolbox/+otp/+oregonator/+presets/EnrightHull.m
+++ b/toolbox/+otp/+oregonator/+presets/EnrightHull.m
@@ -1,25 +1,17 @@
 classdef EnrightHull < otp.oregonator.OregonatorProblem
-    % The Oregonator configuration used in problem CHM 9 of :cite:p:`EH76`. It uses time span $t \in [0, 300]$, initial
+    % The Oregonator configuration used in problem CHM 9 of :cite:p:`EH76`. It uses time span $t ∈ [0, 300]$, initial
     % condition $y_0 = [4, 1.1, 4]^T$, and parameters
     %
     % $$
     % f &= 1, \\
-    % q &= 8.375 \cdot 10^{-6}, \\
+    % q &= 8.375 \times 10^{-6}, \\
     % s &= 77.27, \\
     % w &= 0.1610.
     % $$
 
     methods
         function obj = EnrightHull
-            % Create the EnrightHull problem object.
-            %
-            % Parameters
-            % ----------
-            %
-            % Returns
-            % -------
-            % obj : EnrightHull
-            %    The constructed problem.
+            % Create the EnrightHull Oregonator problem object.
 
             tspan = [0, 300];
             y0 = [4; 1.1; 4];
diff --git a/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m b/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
index eeae9466..7bee036c 100644
--- a/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
+++ b/toolbox/+otp/+oregonator/+presets/GottwaldWanner.m
@@ -1,10 +1,10 @@
 classdef GottwaldWanner < otp.oregonator.OregonatorProblem
-    % The Oregonator configuration from :cite:p:`GW82` which uses time span $t \in [0, 302.85805]$, initial condition
+    % The Oregonator configuration from :cite:p:`GW82` which uses time span $t ∈ [0, 302.85805]$, initial condition
     % $y_0 = [4, 1.331391, 2.852348]^T$, and parameters
     %
     % $$
     % f &= 1, \\
-    % q &= 8.375 \cdot 10^{-6}, \\
+    % q &= 8.375 \times 10^{-6}, \\
     % s &= 77.27, \\
     % w &= 0.1610.
     % $$
@@ -13,15 +13,7 @@
 
     methods
         function obj = GottwaldWanner
-            % Create the GottwaldWanner problem object.
-            %
-            % Parameters
-            % ----------
-            %
-            % Returns
-            % -------
-            % obj : GottwaldWanner
-            %    The constructed problem.
+            % Create the GottwaldWanner Oregonator problem object.
 
             tspan = [0, 302.85805];
             y0 = [4; 1.331391; 2.852348];
diff --git a/toolbox/+otp/+oregonator/OregonatorProblem.m b/toolbox/+otp/+oregonator/OregonatorProblem.m
index d4c57e10..e2a6e49d 100644
--- a/toolbox/+otp/+oregonator/OregonatorProblem.m
+++ b/toolbox/+otp/+oregonator/OregonatorProblem.m
@@ -55,11 +55,7 @@
             %    The initial conditions.
             % parameters : OregonatorParameters
             %    The parameters.
-            %
-            % Returns
-            % -------
-            % obj : OregonatorProblem
-            %    The constructed problem.
+            
             obj@otp.Problem('Oregonator', 3, timeSpan, y0, parameters);
         end
     end
diff --git a/toolbox/+otp/+protherorobinson/+presets/Canonical.m b/toolbox/+otp/+protherorobinson/+presets/Canonical.m
index 5e720645..0171bc71 100644
--- a/toolbox/+otp/+protherorobinson/+presets/Canonical.m
+++ b/toolbox/+otp/+protherorobinson/+presets/Canonical.m
@@ -1,6 +1,6 @@
 classdef Canonical < otp.protherorobinson.ProtheroRobinsonProblem
     % The Prothero–Robinson configuration from :cite:p:`PR74` (p. 159) based on :cite:p:`SLH70` (p. 272). It uses time
-    % span $t \in [0, 15]$, initial condition $y(0) = φ(0)$, and parameters
+    % span $t ∈ [0, 15]$, initial condition $y(0) = φ(0)$, and parameters
     %
     % $$
     % λ &= -200, \\
@@ -19,19 +19,12 @@
             %    - ``Lambda`` – The stiffness parameter and eigenvalue of the Jacobian $λ$.
             %    - ``Phi`` – The exact solution.
             %    - ``DPhi`` – The time derivative of $φ(t)$.
-            %
-            % Returns
-            % -------
-            % obj : ProtheroRobinsonProblem
-            %    The constructed problem
             
             params = otp.protherorobinson.ProtheroRobinsonParameters('Lambda', -200, ...
                 'Phi', @(t) 10 - (10 + t) .* exp(-t), 'DPhi', @(t) (9 + t) .* exp(-t), varargin{:});
             tspan = [0, 15];
             y0 = params.Phi(tspan(1));
-            
             obj = obj@otp.protherorobinson.ProtheroRobinsonProblem(tspan, y0, params);
         end
-        
     end
 end
\ No newline at end of file
diff --git a/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m b/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m
index f3c65f7a..298445e0 100644
--- a/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m
+++ b/toolbox/+otp/+protherorobinson/ProtheroRobinsonProblem.m
@@ -37,11 +37,7 @@
             %    The initial conditions.
             % parameters : ProtheroRobinsonParameters
             %    The parameters.
-            %
-            % Returns
-            % -------
-            % obj : ProtheroRobinsonProblem
-            %    The constructed problem.
+            
             obj@otp.Problem('Prothero–Robinson', [], timeSpan, y0, parameters);
         end
     end
diff --git a/toolbox/+otp/+robertson/+presets/Canonical.m b/toolbox/+otp/+robertson/+presets/Canonical.m
index d769a73c..fbc97da5 100644
--- a/toolbox/+otp/+robertson/+presets/Canonical.m
+++ b/toolbox/+otp/+robertson/+presets/Canonical.m
@@ -1,11 +1,12 @@
 classdef Canonical < otp.robertson.RobertsonProblem
-    % The original configuration :cite:p:`Rob67` with $t \in [0, 40]$, $y_0 = [1, 0, 0]^T$, and parameters
+    % The original configuration :cite:p:`Rob66` with $t ∈ [0, 40]$, $y_0 = [1, 0, 0]^T$, and parameters
     %
     % $$
-    % K_1 &= 4\times 10^{-2}, \\
-    % K_2 &= 3\times 10^7, \\
-    % K_3 &= 10^4.
+    % k_1 &= 4 \times 10^{-2}, \\
+    % k_2 &= 3 \times 10^7, \\
+    % k_3 &= 10^4.
     % $$
+
     methods
         function obj = Canonical(varargin)
             % Create the Canonical Robertson problem object.
@@ -15,21 +16,14 @@
             % varargin
             %    A variable number of name-value pairs. The accepted names are
             %
-            %    - ``K1`` – Value of $K_1$.
-            %    - ``K2`` – Value of $K_2$.
-            %    - ``K3`` – Value of $K_3$.
-            %
-            % Returns
-            % -------
-            % obj : RobertsonProblem
-            %    The constructed problem.
+            %    - ``K1`` – Value of $k_1$.
+            %    - ``K2`` – Value of $k_2$.
+            %    - ``K3`` – Value of $k_3$.
             
             y0 = [1; 0; 0];
             tspan = [0, 40];
             params = otp.robertson.RobertsonParameters('K1', 0.04, 'K2', 3e7, 'K3', 1e4, varargin{:});
-            
             obj = obj@otp.robertson.RobertsonProblem(tspan, y0, params);
         end
-
     end
 end
diff --git a/toolbox/+otp/+robertson/+presets/ROBER.m b/toolbox/+otp/+robertson/+presets/ROBER.m
index 3df1ed30..a77d054a 100644
--- a/toolbox/+otp/+robertson/+presets/ROBER.m
+++ b/toolbox/+otp/+robertson/+presets/ROBER.m
@@ -2,21 +2,15 @@
     % The Robertson configuration from :cite:p:`HW96` (p. 144). This differs from :class:`Canonical` only in the time
     % span which is extended to $[0, 10^{11}]$. This presents a challenge for numerical integrators to preserve solution
     % positivity.
+
     methods
         function obj = ROBER
             % Create the ROBER Robertson problem object.
-            %
-            % Returns
-            % -------
-            % obj : RobertsonProblem
-            %    The constructed problem.
             
             y0 = [1; 0; 0];
-            tspan = [0; 1e11];
+            tspan = [0, 1e11];
             params = otp.robertson.RobertsonParameters('K1', 0.04, 'K2', 3e7, 'K3', 1e4);
-            
             obj = obj@otp.robertson.RobertsonProblem(tspan, y0, params);
         end
-
     end
 end
diff --git a/toolbox/+otp/+robertson/RobertsonProblem.m b/toolbox/+otp/+robertson/RobertsonProblem.m
index 64a7e186..8564a297 100644
--- a/toolbox/+otp/+robertson/RobertsonProblem.m
+++ b/toolbox/+otp/+robertson/RobertsonProblem.m
@@ -1,28 +1,27 @@
 classdef RobertsonProblem < otp.Problem
     % A simple, stiff chemical reaction.
     % 
-    % The Robertson problem :cite:p:`Rob67` is a simple, stiff chemical reaction that
-    % models the concentrations of chemical species $A$, $B$, and $C$ 
-    % involved in the reactions
+    % The Robertson problem :cite:p:`Rob66` is a simple, stiff chemical reaction that models the concentrations of
+    % chemical species $\ce{A}$, $\ce{B}$, and $\ce{C}$ involved in the reactions
+    %
     % $$
-    % \begin{align*}
-    % &A  \rightarrow B      &~~\text{at rate}~~ &K_1, \\
-    % &B + B \rightarrow C + B  &~~\text{at rate}~~  &K_2, \\
-    % &B + C \rightarrow A + C  &~~\text{at rate}~~  &K_3.
-    % \end{align*}
+    % \ce{
+    % A &->[$k_1$] B \\
+    % B + B &->[$k_2$] C + B \\
+    % B + C &->[$k_3$] A + C
+    % }
     % $$
+    %
     % These correspond to the ODE system, 
     % 
     % $$ 
-    % \begin{align*}
-    % y_1' &= -K_1 y_1 + K_3  y_2 y_3, \\
-    % y_2' &= K_1 y_1 - K_2 y_2^2 - K_3 y_2 y_3, \\
-    % y_3' &= K_2 y_2^2.
-    % \end{align*}
+    % y_1' &= -k_1 y_1 + k_3  y_2 y_3, \\
+    % y_2' &= k_1 y_1 - k_2 y_2^2 - k_3 y_2 y_3, \\
+    % y_3' &= k_2 y_2^2.
     % $$
-    % The reaction rates $K_1$, $K_2$, and $K_3$ often range from slow to very fast
-    % which makes the problem challenging. This has made it a popular test for
-    % implicit time-stepping schemes.
+    %
+    % The reaction rates $k_1$, $k_2$, and $k_3$ often range from slow to very fast which makes the problem challenging.
+    % This has made it a popular test for implicit time-stepping schemes.
     % 
     % Notes
     % -----
@@ -31,12 +30,11 @@
     % +---------------------+-------------------------------------------------+
     % | Number of Variables | 3                                               |
     % +---------------------+-------------------------------------------------+
-    % | Stiff               | typically, depending on $K_1$, $K_2$, and $K_3$ |
+    % | Stiff               | typically, depending on $k_1$, $k_2$, and $k_3$ |
     % +---------------------+-------------------------------------------------+
     %
     % Example
     % -------
-    %
     % >>> problem = otp.robertson.presets.Canonical;
     % >>> problem.Parameters.K1 = 100;
     % >>> problem.Parameters.K2 = 1000;
@@ -60,11 +58,7 @@
             %    The initial conditions.
             % parameters : RobertsonParameters
             %    The parameters.
-            %
-            % Returns
-            % -------
-            % obj : RobertsonProblem
-            %    The constructed problem.
+            
             obj@otp.Problem('Robertson', 3, timeSpan, y0, parameters);
         end
     end
diff --git a/toolbox/+otp/RHS.m b/toolbox/+otp/RHS.m
index 869bf940..32324d92 100644
--- a/toolbox/+otp/RHS.m
+++ b/toolbox/+otp/RHS.m
@@ -30,6 +30,8 @@
     %
     % See Also
     % --------
+    % otp.Problem.RHS
+    %
     % odeset : https://www.mathworks.com/help/matlab/ref/odeset.html
 
     properties (SetAccess = private)

From 4a08b8c5c216b5e937c08fd029374ee9c688c6da Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Sun, 5 Nov 2023 22:07:49 -0800
Subject: [PATCH 17/20] Fix capitalization

---
 toolbox/+otp/+robertson/RobertsonParameters.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/toolbox/+otp/+robertson/RobertsonParameters.m b/toolbox/+otp/+robertson/RobertsonParameters.m
index fea738ea..6ed8e2cc 100644
--- a/toolbox/+otp/+robertson/RobertsonParameters.m
+++ b/toolbox/+otp/+robertson/RobertsonParameters.m
@@ -1,11 +1,11 @@
 classdef RobertsonParameters < otp.Parameters
     % Parameters for the Robertson problem.
     properties
-        % The reaction rate $K_1$.
+        % The reaction rate $k_1$.
         K1 %MATLAB ONLY: (1, 1) {mustBeReal, mustBeNonnegative}
-        % The reaction rate $K_2$.
+        % The reaction rate $k_2$.
         K2 %MATLAB ONLY: (1, 1) {mustBeReal, mustBeNonnegative}
-        % The reaction rate $K_3$.
+        % The reaction rate $k_3$.
         K3 %MATLAB ONLY: (1, 1) {mustBeReal, mustBeNonnegative}
     end
 

From a2a3fa3ecdcfc54362b54b4f0fe533eedce71140 Mon Sep 17 00:00:00 2001
From: Andrey A Popov <54905478+AndreyAPopov@users.noreply.github.com>
Date: Tue, 7 Nov 2023 22:07:44 -0600
Subject: [PATCH 18/20] fixes for the project file (#123)

---
 toolboxPackaging.prj | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/toolboxPackaging.prj b/toolboxPackaging.prj
index 6d548e06..355c9134 100644
--- a/toolboxPackaging.prj
+++ b/toolboxPackaging.prj
@@ -7,7 +7,7 @@
     <param.summary>A MATLAB suite of initial value problems</param.summary>
     <param.description />
     <param.screenshot>${PROJECT_ROOT}/images/logo.png</param.screenshot>
-    <param.version>0.0.1</param.version>
+    <param.version>0.2.0</param.version>
     <param.output>${PROJECT_ROOT}/ODE Test Problems.mltbx</param.output>
     <param.products.name />
     <param.products.id />
@@ -44,7 +44,6 @@
     <unset>
       <param.email />
       <param.description />
-      <param.screenshot />
       <param.output />
       <param.products.name />
       <param.products.id />

From 29e48590aa05c6f01409ab719bd6995e9b1c47e1 Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Tue, 7 Nov 2023 20:08:48 -0800
Subject: [PATCH 19/20] Bumb version

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 087450f2..0879794c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,5 +1,5 @@
 Name: ODE Test Problems
-Version: 0.0.1
+Version: 0.2.0
 Date: May 10, 2023
 depends: octave (>= 6.4.0)
 Author: Steven Roberts, Andrey A. Popov, Arash Sarshar, Adrian Sandu

From b040d626231ec7d4d54719b4a2732161e6ac29bb Mon Sep 17 00:00:00 2001
From: Steven Roberts <sroberts994@gmail.com>
Date: Tue, 7 Nov 2023 20:09:27 -0800
Subject: [PATCH 20/20] Update DESCRIPTION

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 0879794c..44a824bd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Name: ODE Test Problems
 Version: 0.2.0
-Date: May 10, 2023
+Date: November 7, 2023
 depends: octave (>= 6.4.0)
 Author: Steven Roberts, Andrey A. Popov, Arash Sarshar, Adrian Sandu
 Maintainer: Steven Roberts