diff --git a/docs/references.bib b/docs/references.bib
index 728b7b8a..a8c991d8 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -178,4 +178,22 @@ @article{Pet86
  year = {1986},
  url = {https://www.jstor.org/stable/2157625}
 }
-
+@techreport{LS98,
+  title        = {Test Set for Initial Value Problem Solvers},
+  author       = {W.M. Lioen and J.J.B. de Swart},
+  institution  = {Modelling, Analysis and Simulation (MAS)},
+  number       = {MAS-R9832},
+  year         = {1998},
+  month        = {dec},
+  type         = {Report},
+  address      = {CWI, P.O. Box 94079, 1090 GB Amsterdam, The Netherlands},
+  issn         = {1386-3703}
+}
+@phdthesis{Sto98,
+  author       = {Walter Johannes Henricus Stortelder},
+  title        = {Parameter Estimation in Nonlinear Dynamical Systems},
+  year         = {1998},
+  school       = {Centrum Wiskunde \& Informatica},
+  address      = {Amsterdam, The Netherlands},
+  type         = {PhD thesis}
+}
\ No newline at end of file
diff --git a/toolbox/+otp/+zlakinetics/+presets/Canonical.m b/toolbox/+otp/+zlakinetics/+presets/Canonical.m
index 3b7b4973..bd33e300 100644
--- a/toolbox/+otp/+zlakinetics/+presets/Canonical.m
+++ b/toolbox/+otp/+zlakinetics/+presets/Canonical.m
@@ -1,11 +1,8 @@
-classdef Canonical < otp.zlakinetics.ZLAKineticsProblem
-    %CANONICAL The problem as described in the literature
-    %
-    % See:
-    %    Walter Johannes Henricus Stortelder. Parameter estimation in nonlinear dynamical
-    %    systems. PhD thesis, Centrum Wiskunde & Informatica, 
-    %    Amsterdam, The Netherlands, 1998.
-    %
+classdef Canonical < otp.zlakinetics.ZLAKineticsProblem 
+    % Canonical preset for the ZLA kinetics problem as descrobe in :cite:p:`Sto98` (pg. 118)
+    % where the initial conditions are $y_0 = [0.444, 0.00123, 0, 0.007, 0, 0]^T$ is also
+    % provided. The time span is $[0, 180]$.
+
     methods
         function obj = Canonical(varargin)
             tspan = [0, 180];
diff --git a/toolbox/+otp/+zlakinetics/ZLAKineticsParameters.m b/toolbox/+otp/+zlakinetics/ZLAKineticsParameters.m
index 3777f30f..eab48d59 100644
--- a/toolbox/+otp/+zlakinetics/ZLAKineticsParameters.m
+++ b/toolbox/+otp/+zlakinetics/ZLAKineticsParameters.m
@@ -1,23 +1,25 @@
 classdef ZLAKineticsParameters < otp.Parameters
     % Parameters for the ZLA kinetics problem.
+    %
+
     
     properties
-        %k are the reaction rates of ancillary reactions
+        % The reaction rates of ancillary reactions.
         k %MATLAB ONLY: (1,4) {otp.utils.validation.mustBeNumerical} = [18.7, 0.58, 0.09, 0.42];
 
-        %K is the ZLA reaction rate
+        %K The ZLA reaction rate.
         K %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} = 34.4;
 
-        %KLA is the mass transfer coefficient
+        % The mass transfer coefficient.
         KlA %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} = 3.3;
 
-        %KS is the equilibrium constant
+        % The equilibrium constant.
         Ks %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} = 115.83;
 
-        %PCO2 is the partial CO_2 pressure
+        % Partial $\text{CO}_2$ pressure.
         PCO2 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} = 0.9;
 
-        %H is the Henry constant
+        % The Henry constant.
         H %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical} = 737;
     end
 
diff --git a/toolbox/+otp/+zlakinetics/ZLAKineticsProblem.m b/toolbox/+otp/+zlakinetics/ZLAKineticsProblem.m
index f6d22a60..df65d489 100644
--- a/toolbox/+otp/+zlakinetics/ZLAKineticsProblem.m
+++ b/toolbox/+otp/+zlakinetics/ZLAKineticsProblem.m
@@ -1,4 +1,52 @@
 classdef ZLAKineticsProblem < otp.Problem
+    % ZLA Kinetics model also known as The Chemical Akzo Nobel problem is a differential
+    % equation that describes the interaction of fictitious species, FLB and ZHU to form two other 
+    % species, ZLA and ZLH. The complete reaction is given in :cite:p:`Sto98`.
+    %
+    % Modeling the concentrations of the species involved in the reaction as $y_{i=\{1,\ldots, 6 \}}$ 
+    % and various reaction rates as $r_{i=\{1,\ldots, 5 \}}$ leads to the following system of ODEs:
+    %
+    % $$
+    % \begin{aligned}
+    % & y_1^{\prime}=-2 r_1 + r_2 - r_3 - r_4, \\
+    % & y_2^{\prime}=-\frac{1}{2} r_1 - r_4 - \frac{1}{2} r_5 + F_{\text{in}}, \\
+    % & y_3^{\prime}=r_1 - r_2 + r_3, \\
+    % & y_4^{\prime}=-r_2 + r_3 - 2 r_4, \\
+    % & y_5^{\prime}=r_2 - r_3 + r_5, \\
+    % & y_6^{\prime}=K_s \cdot y_1 \cdot y_2 \cdot y_6,\\
+    % & F_{\text {in }} = \text{klA} \cdot\left(\frac{\text{pCO}_2}{H}-y_2\right).
+    % \end{aligned}
+    % $$
+    % and the reaction rates are given by:
+    % $$
+    % \begin{aligned}
+    % r_1 & =k_1 \cdot y_1^4 \cdot y_2^{\frac{1}{2}}, \\
+    % r_2 & =k_2 \cdot y_3 \cdot y_4, \\
+    % r_3 & =\frac{k_2}{K} \cdot y_1 \cdot y_5, \\
+    % r_4 & =k_3 \cdot y_1 \cdot y_4^2, \\
+    % r_5 & =k_4 \cdot y_6^2 \cdot y_2^{\frac{1}{2}}.
+    % \end{aligned}
+    % $$
+    % 
+    % The parameters $k = [k_i]_{i=\{1,\ldots, 4 \}}$, $K$, $\text{klA}$, $K_s$, $\text{pCO}_2$, and $H$ are the rates and constants of the reaction.
+    %
+    % Notes
+    % -----
+    % +---------------------+-----------------------------------------+
+    % | Type                | ODE                                     |
+    % +---------------------+-----------------------------------------+
+    % | Number of Variables | 6                                       |
+    % +---------------------+-----------------------------------------+
+    % | Stiff               | yes, depending on parameters            |
+    % +---------------------+-----------------------------------------+
+    %
+    % Example
+    % -------
+    % >>> problem = otp.zlakinetics.presets.Canonical;
+    % >>> problem.TimeSpan(end) = 15;
+    % >>> sol = problem.solve();
+    % >>> problem.plot(sol);
+    
     methods
         function obj = ZLAKineticsProblem(timeSpan, y0, parameters)
             obj@otp.Problem('ZLA-Kinetics', 6, timeSpan, y0, parameters);