From 7bf8f805dc4d60260310d0d058eab4009a4d3c47 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 11 Jul 2024 14:17:54 +0200
Subject: [PATCH 01/17] Updated basic examples

---
 docs/src/index.md | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 50c21d4e..22d89feb 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -78,8 +78,8 @@ To solve this PDS together with initial values ``u_1(0)=u_2(0)=2`` on the time d
 ```@example LotkaVolterra
 using PositiveIntegrators # load PDSProblem
 
-P(u, p, t) = [2*u[1]  0.0; u[1]*u[2]  0.0] # Production matrix
-d(u, p, t) = [0.0; u[2]] # Destruction vector
+P(u, p, t) = [2*u[1]  0; u[1]*u[2]  0] # Production matrix
+d(u, p, t) = [0; u[2]] # Destruction vector
 
 u0 = [2.0; 2.0] # initial values
 tspan = (0.0, 10.0) # time span
@@ -88,12 +88,10 @@ tspan = (0.0, 10.0) # time span
 prob = PDSProblem(P, d, u0, tspan)
 nothing #hide
 ```
-Now that the problem has been created, we can solve it with any of the methods of [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Here we use the method `Tsit5()`. Please note that [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) currently only provides methods for positive and conservative PDS, see below.
+Now that the problem has been created, we can solve it with any method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl). In the following, we use the method `MPRK22(1.0)`. In addition, we could also use any method provided by [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/), but these might possibly generate negative approximations.
 
 ```@example LotkaVolterra
-using OrdinaryDiffEq  #load Tsit5
-
-sol = solve(prob, Tsit5())
+sol = solve(prob, MPRK22(1.0))
 nothing # hide
 ```
 Finally, we can use [Plots.jl](https://docs.juliaplots.org/stable/) to visualize the solution.
@@ -119,7 +117,6 @@ This shows that the sum of the state variables of a conservative PDS remains con
 \sum_{i=1}^N y_i(t) = \sum_{i=1}^N y_i(0)
 ```
 for all times ``t>0``.
-Moreover, a conservative PDS is completely defined by the square matrix ``\mathbf P=(p_{ij})_{i,j=1,\dots,N}``. There is no need to store an additional vector of destruction terms since ``d_{ij} = p_{ji}`` for all ``i,j=1,\dots,N``.
 
 One specific example of a conservative PDS is the SIR model
 ```math
@@ -165,7 +162,7 @@ tspan = (0.0, 100.0); # time span
 prob = ConservativePDSProblem(P, u0, tspan)
 nothing # hide
 ```
-Since the SIR model is not only conservative but also positive, we can use any MPRK scheme from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) to solve it. Here we use `MPRK22(1.0)`.
+Since the SIR model is not only conservative but also positive, we can use any scheme from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) to solve it. Here we use `MPRK22(1.0)`.
 Please note that any method from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) can be used as well, but might possibly generate negative approximations.
 
 ```@example SIR
@@ -176,9 +173,10 @@ Finally, we can use [Plots.jl](https://docs.juliaplots.org/stable/) to visualize
 ```@example SIR
 using Plots
 
-plot(sol, legend=:right)
+plot(sol, label = ["S" "I" "R"], legend=:right)
+plot!(sol, idxs = ((t, S, I, R) -> (t, S + I + R), 0, 1, 2, 3), label = "S+I+R") #Plot S+I+R over time.
 ```
-
+We see that there is always a nonnegative number of people in each compartment, while the population ``S+I+R`` remains constant over time.
 
 ## Referencing
 

From 4009c33e927b5d394c2ab59d073d0fc50535cb0b Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Thu, 11 Jul 2024 16:47:34 +0200
Subject: [PATCH 02/17] Edited doc strings for MPRK and SSPMPRK schemes

---
 src/mprk.jl    | 35 ++++++++++++++++++++++++-----------
 src/sspmprk.jl | 10 ++++++++--
 2 files changed, 32 insertions(+), 13 deletions(-)

diff --git a/src/mprk.jl b/src/mprk.jl
index 8b327b0b..9434cec6 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -122,13 +122,14 @@ The first-order modified Patankar-Euler algorithm for production-destruction sys
 first-order accurate, unconditionally positivity-preserving, and
 linearly implicit.
 
-The scheme was introduced by Burchard et al for conservative production-destruction systems.
+The scheme was introduced by Burchard et al. for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 
-``u_i^{n+1} = u_i^n + Δt \\sum_{j, j≠i} \\biggl(p_{ij}^n \\frac{u_j^{n+1}}{u_j^n}-d_{ij}^n \\frac{u_i^{n+1}}{u_i^n}\\biggr) + {\\Delta}t p_{ii}^n - Δt d_{ii}^n\\frac{u_i^{n+1}}{u_i^n}``.
+``u_i^{n+1} = u_i^n + Δt \\sum_{j, j≠i} \\biggl(p_{ij}^n \\frac{u_j^{n+1}}{u_j^n}-d_{ij}^n \\frac{u_i^{n+1}}{u_i^n}\\biggr) + {\\Delta}t p_{ii}^n - Δt d_{ii}^n\\frac{u_i^{n+1}}{u_i^n}``,
 
-The modified Patankar-Euler method requires the special structure of a
-[`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref).
+where ``p_{ij}^n = p_{ij}(t^n,\\mathbf u^n)`` and ``d_{ij}^n = d_{ij}(t^n,\\mathbf u^n)``.
+
+The modified Patankar-Euler method requires the special structure of a [`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref).
 
 You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
@@ -137,6 +138,8 @@ You can also choose the parameter `small_constant` which is added to all Patanka
 to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
 `floatmin` of the floating point type used.
 
+The current implementation only supports fixed time steps. 
+
 ## References
 
 - Hans Burchard, Eric Deleersnijder, and Andreas Meister.
@@ -329,14 +332,18 @@ end
     MPRK22(α; [linsolve = ..., small_constant = ...])
 
 A family of second-order modified Patankar-Runge-Kutta algorithms for
-production-destruction systems. Each member of this family is an one-step, two-stage method which is
+production-destruction systems. Each member of this family is an adaptive, one-step, two-stage method which is
 second-order accurate, unconditionally positivity-preserving, and linearly
-implicit. The parameter `α` is described by Kopecz and Meister (2018) and
+implicit. In this implementation the stage-values are conservative as well. 
+The parameter `α` is described by Kopecz and Meister (2018) and
 studied by Izgin, Kopecz and Meister (2022) as well as
-Torlo, Öffner and Ranocha (2022).
+Torlo, Öffner and Ranocha (2022).  
+
+This method supports adaptive time stepping, using the Patankar-weight denominators 
+``σ_i``, see Kopecz and Meister (2018), as first order approximations to estimate the error.
 
 The scheme was introduced by Kopecz and Meister for conservative production-destruction systems.
-For nonconservative production–destruction systems we use the straight forward extension
+For nonconservative production–destruction systems we use a straight forward extension
 analogous to [`MPE`](@ref).
 
 This modified Patankar-Runge-Kutta method requires the special structure of a
@@ -712,12 +719,15 @@ end
 
 A family of third-order modified Patankar-Runge-Kutta schemes for
 production-destruction systems, which is based on the two-parameter family of third order explicit Runge--Kutta schemes.
-Each member of this family is a one-step method with four-stages which is
+Each member of this family is an adaptive, one-step method with four-stages which is
 third-order accurate, unconditionally positivity-preserving, conservative and linearly
 implicit. In this implementation the stage-values are conservative as well.
 The parameters `α` and `β` must be chosen such that the Runge--Kutta coefficients are nonnegative,
 see Kopecz and Meister (2018) for details.
 
+These methods support adaptive time stepping, using the Patankar-weight denominators 
+``σ_i``, see Kopecz and Meister (2018), as second order approximations to estimate the error.
+
 The scheme was introduced by Kopecz and Meister for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 analogous to [`MPE`](@ref).
@@ -809,15 +819,18 @@ end
 """
     MPRK43II(γ; [linsolve = ..., small_constant = ...])
 
-A family of third-order modified Patankar-Runge-Kutta schemes for (conservative)
+A family of third-order modified Patankar-Runge-Kutta schemes for
 production-destruction systems, which is based on the one-parameter family of third order explicit Runge--Kutta schemes with
 non-negative Runge--Kutta coefficients.
-Each member of this family is a one-step method with four stages which is
+Each member of this family is an adaptive, one-step method with four stages which is
 third-order accurate, unconditionally positivity-preserving, conservative and linearly
 implicit. In this implementation the stage-values are conservative as well. The parameter `γ` must satisfy
 `3/8 ≤ γ ≤ 3/4`.
 Further details are given in Kopecz and Meister (2018).
 
+This method supports adaptive time stepping, using the Patankar-weight denominators 
+``σ_i``, see Kopecz and Meister (2018), as second order approximations to estimate the error.
+
 The scheme was introduced by Kopecz and Meister for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 analogous to [`MPE`](@ref).
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index cb2f1144..bd871d33 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -3,7 +3,7 @@
     SSPMPRK22(α, β; [linsolve = ..., small_constant = ...])
 
 A family of second-order modified Patankar-Runge-Kutta algorithms for
-production-destruction systems. Each member of this family is a one-step, two-stage method which is
+production-destruction systems. Each member of this family is an adaptive, one-step, two-stage method which is
 second-order accurate, unconditionally positivity-preserving, and linearly
 implicit. The parameters `α` and `β` are described by Huang and Shu (2019) and
 studied by Huang, Izgin, Kopecz, Meister and Shu (2023).
@@ -11,6 +11,10 @@ The difference to [`MPRK22`](@ref) is that this method is based on the SSP formu
 an explicit second-order Runge-Kutta method. This family of schemes contains the [`MPRK22`](@ref) family,
 where `MPRK22(α) = SSMPRK22(0, α)` applies.
 
+This method supports adaptive time stepping, using the first order approximations 
+``(σ_i - u_i^n) / τ + u_i^n`` with ``τ=1+(α_{21}β_{10}^2)/(β_{20}+β_{21})``, 
+see (2.7) in Huang and Shu (2019), to estimate the error.
+
 The scheme was introduced by Huang and Shu for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 analogous to [`MPE`](@ref).
@@ -408,7 +412,7 @@ end
     SSPMPRK43([linsolve = ..., small_constant = ...])
 
 A third-order modified Patankar-Runge-Kutta algorithm for
-production-destruction systems. This scheme is a one-step, two-stage method which is
+production-destruction systems. This scheme is a one-step, four-stage method which is
 third-order accurate, unconditionally positivity-preserving, and linearly
 implicit. The scheme is described by Huang, Zhao and Shu (2019) and
 studied by Huang, Izgin, Kopecz, Meister and Shu (2023).
@@ -429,6 +433,8 @@ You can also choose the parameter `small_constant` which is added to all Patanka
 to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
 `floatmin` of the floating point type used.
 
+The current implementation only supports fixed time steps. 
+
 ## References
 
 - Juntao Huang, Weifeng Zhao and Chi-Wang Shu.

From 73bfa3aa63a6597c18be41ad3fc949c85a760359 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 12 Jul 2024 16:42:23 +0200
Subject: [PATCH 03/17] Started tutorial for linear advection

---
 docs/make.jl                 |  1 +
 docs/src/linear_advection.md | 96 ++++++++++++++++++++++++++++++++++++
 2 files changed, 97 insertions(+)
 create mode 100644 docs/src/linear_advection.md

diff --git a/docs/make.jl b/docs/make.jl
index 56c5d6a6..467a19d8 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -73,6 +73,7 @@ makedocs(modules = [PositiveIntegrators],
          # Explicitly specify documentation structure
          pages = [
              "Home" => "index.md",
+             "Tutorials" => "linear_advection.md",
              "API reference" => "api_reference.md",
              "Contributing" => "contributing.md",
              "Code of conduct" => "code_of_conduct.md",
diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
new file mode 100644
index 00000000..8e5b65a7
--- /dev/null
+++ b/docs/src/linear_advection.md
@@ -0,0 +1,96 @@
+# Tutorial: Solution of the linear advection equation
+
+This tutorial is about the efficient solution of production-destruction systems (PDS) with a large number of differential equations. 
+We will explore several ways to represent such large systems and assess their efficiency. 
+
+## Definition of the production-destruction system
+
+One example of the occurrence of a PDS with a large number of equations is the space discretization of a partial differential equation. In this tutorial we want to solve the linear advection equation
+
+$$\partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x)$$
+
+with $a>0$, $t≥ 0$, $x\in[0,1]$ and periodic boundary condtions. To keep things as simple as possible, we 
+discretize the space domain as $0=x_0<x_1\dots <x_{N-1}<x_N=1$ with $x_i = i Δ x$ for $i=0,\dots,N$ and $Δx=1/N$. An upwind discretization of the spatial derivative yields the ODE system
+
+$$\begin{aligned}
+&\partial_t u_1(t) =-\frac{a}{Δx}\bigl(u_1(t)-u_{N}(t)\bigr),\\
+&\partial_t u_i(t) =-\frac{a}{Δx}\bigl(u_i(t)-u_{i-1}(t)\bigr),\quad i=2,\dots,N,
+\end{aligned}
+$$ 
+
+where $u_i(t)$ is an approximation of $u(t,x_i)$ for $i=1,\dots, N$.
+This system can also be written as $\partial_t \mathbf u(t)=\mathbf A\mathbf u(t)$ with $\mathbf u(t)=(u_1(t),\dots,u_N(t))$ and 
+
+$$
+\mathbf A= \frac{a}{Δ x}\begin{bmatrix}-1&0&\dots&0&1\\1&-1&\ddots&&0\\0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\0&\dots&0&1&-1\end{bmatrix}.
+$$
+
+In particular the matrix $\mathbf A$ shows, that there is a single production term and a single destruction term per equation. 
+Furthermore, the system is conservative as $\mathbf A$ has column sum zero.
+To be precise, the production matrix $\mathbf P = (p_{i,j})$ of this conservative PDS is given by
+
+$$\begin{aligned}
+&p_{1,N}(t,\mathbf u(t)) = \frac{a}{Δ x}u_N(t),\\
+&p_{i,i-1}(t,\mathbf u(t)) = \frac{a}{Δ x}u_{i-1}(t),\quad i=2,\dots,N.
+\end{aligned}
+$$
+
+Since the PDS is conservative, we have $d_{i,j}=p_{j,i}$ and the system is fully determined by the production matrix $\mathbf P$.
+
+## Solution of the production-destruction system
+
+Now we are ready to define a `ConservativePDSProblem` and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). In the following we use $a=1$, $N=1000$ and the time domain $t\in[0,1]$. Moreover, we choose the step function
+
+$$
+u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
+$$
+
+as initial condition. Due to the peridoic boundary conditions and the transport velocity $a=1$, the solution at time $t=1$ is identical to the inital distribution, i.e. $u(1,x) = u_0(x)$.
+
+```@setup LinearAdvection
+import Pkg; Pkg.add("PositiveIntegrators"); Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots")
+```
+```@example LinearAdvection
+N = 1000 # number of subintervals
+dx = 1/N; # mesh width
+x = LinRange(dx, 1.0, N+1) # discretization points x_1,...,x_N = x_0
+u0 = 0.0 .+ (0.4 .≤ x .≤ 0.6) .* 1.0 # initial solution
+tspan = (0.0, 1.0) # time domain
+
+nothing #hide
+```
+
+As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are
+1. an in-place implementation with a dense matrix,
+2. an in-place implementation with a sparse matrix.
+
+### Standard in-place implementation
+
+```@example LinearAdvection
+using PositiveIntegrators # load ConservativePDSProblem
+
+function lin_adv_P!(P, u, p, t)
+    P .= 0.0
+    N = length(u)
+    dx = 1 / N
+    P[1, N] = u[N] / dx
+    for i in 2:N
+        P[i, i - 1] = u[i - 1] / dx
+    end
+    return nothing
+end
+
+prob = ConservativePDSProblem(lin_adv_P!, u0, tspan) # create the PDS
+
+sol = solve(prob, MPRK43I(1.0, 0.5); save_everystep = false)
+
+nothing #hide
+```
+
+```@example LinearAdvection
+using Plots
+
+plot(x,u0)
+plot!(x, last(sol.u))
+```
+

From 6769b68cf1c4b557e8fc687530031d73d3cea1fe Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 12 Jul 2024 16:46:28 +0200
Subject: [PATCH 04/17] spell check

---
 docs/src/linear_advection.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 8e5b65a7..9648efcf 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -9,7 +9,7 @@ One example of the occurrence of a PDS with a large number of equations is the s
 
 $$\partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x)$$
 
-with $a>0$, $t≥ 0$, $x\in[0,1]$ and periodic boundary condtions. To keep things as simple as possible, we 
+with $a>0$, $t≥ 0$, $x\in[0,1]$ and periodic boundary conditions. To keep things as simple as possible, we 
 discretize the space domain as $0=x_0<x_1\dots <x_{N-1}<x_N=1$ with $x_i = i Δ x$ for $i=0,\dots,N$ and $Δx=1/N$. An upwind discretization of the spatial derivative yields the ODE system
 
 $$\begin{aligned}
@@ -45,7 +45,7 @@ $$
 u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
 $$
 
-as initial condition. Due to the peridoic boundary conditions and the transport velocity $a=1$, the solution at time $t=1$ is identical to the inital distribution, i.e. $u(1,x) = u_0(x)$.
+as initial condition. Due to the periodic boundary conditions and the transport velocity $a=1$, the solution at time $t=1$ is identical to the initial distribution, i.e. $u(1,x) = u_0(x)$.
 
 ```@setup LinearAdvection
 import Pkg; Pkg.add("PositiveIntegrators"); Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots")

From 68615ad9c826c81503904f0a9a46eda0c6e0d2f6 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 12 Jul 2024 17:00:19 +0200
Subject: [PATCH 05/17] documenter latex

---
 docs/src/linear_advection.md | 42 ++++++++++++++++++++----------------
 1 file changed, 23 insertions(+), 19 deletions(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 9648efcf..9983b48b 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -7,45 +7,49 @@ We will explore several ways to represent such large systems and assess their ef
 
 One example of the occurrence of a PDS with a large number of equations is the space discretization of a partial differential equation. In this tutorial we want to solve the linear advection equation
 
-$$\partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x)$$
+```math
+\partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x)
+```
 
-with $a>0$, $t≥ 0$, $x\in[0,1]$ and periodic boundary conditions. To keep things as simple as possible, we 
-discretize the space domain as $0=x_0<x_1\dots <x_{N-1}<x_N=1$ with $x_i = i Δ x$ for $i=0,\dots,N$ and $Δx=1/N$. An upwind discretization of the spatial derivative yields the ODE system
+with ``a>0``, ``t≥ 0``, ``x\in[0,1]`` and periodic boundary conditions. To keep things as simple as possible, we 
+discretize the space domain as ``0=x_0<x_1\dots <x_{N-1}<x_N=1`` with ``x_i = i Δ x`` for ``i=0,\dots,N`` and ``Δx=1/N``. An upwind discretization of the spatial derivative yields the ODE system
 
-$$\begin{aligned}
+```math
+\begin{aligned}
 &\partial_t u_1(t) =-\frac{a}{Δx}\bigl(u_1(t)-u_{N}(t)\bigr),\\
 &\partial_t u_i(t) =-\frac{a}{Δx}\bigl(u_i(t)-u_{i-1}(t)\bigr),\quad i=2,\dots,N,
 \end{aligned}
-$$ 
+```
 
-where $u_i(t)$ is an approximation of $u(t,x_i)$ for $i=1,\dots, N$.
-This system can also be written as $\partial_t \mathbf u(t)=\mathbf A\mathbf u(t)$ with $\mathbf u(t)=(u_1(t),\dots,u_N(t))$ and 
+where ``u_i(t)`` is an approximation of ``u(t,x_i)`` for ``i=1,\dots, N``.
+This system can also be written as ``\partial_t \mathbf u(t)=\mathbf A\mathbf u(t)`` with ``\mathbf u(t)=(u_1(t),\dots,u_N(t))`` and 
 
-$$
+```math
 \mathbf A= \frac{a}{Δ x}\begin{bmatrix}-1&0&\dots&0&1\\1&-1&\ddots&&0\\0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\0&\dots&0&1&-1\end{bmatrix}.
-$$
+```
 
-In particular the matrix $\mathbf A$ shows, that there is a single production term and a single destruction term per equation. 
-Furthermore, the system is conservative as $\mathbf A$ has column sum zero.
-To be precise, the production matrix $\mathbf P = (p_{i,j})$ of this conservative PDS is given by
+In particular the matrix ``\mathbf A`` shows, that there is a single production term and a single destruction term per equation. 
+Furthermore, the system is conservative as ``\mathbf A`` has column sum zero.
+To be precise, the production matrix ``\mathbf P = (p_{i,j})`` of this conservative PDS is given by
 
-$$\begin{aligned}
+```math
+\begin{aligned}
 &p_{1,N}(t,\mathbf u(t)) = \frac{a}{Δ x}u_N(t),\\
 &p_{i,i-1}(t,\mathbf u(t)) = \frac{a}{Δ x}u_{i-1}(t),\quad i=2,\dots,N.
 \end{aligned}
-$$
+```
 
-Since the PDS is conservative, we have $d_{i,j}=p_{j,i}$ and the system is fully determined by the production matrix $\mathbf P$.
+Since the PDS is conservative, we have ``d_{i,j}=p_{j,i}`` and the system is fully determined by the production matrix ``\mathbf P``.
 
 ## Solution of the production-destruction system
 
-Now we are ready to define a `ConservativePDSProblem` and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). In the following we use $a=1$, $N=1000$ and the time domain $t\in[0,1]$. Moreover, we choose the step function
+Now we are ready to define a `ConservativePDSProblem` and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). In the following we use ``a=1``, ``N=1000`` and the time domain ``t\in[0,1]``. Moreover, we choose the step function
 
-$$
+```math
 u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
-$$
+```
 
-as initial condition. Due to the periodic boundary conditions and the transport velocity $a=1$, the solution at time $t=1$ is identical to the initial distribution, i.e. $u(1,x) = u_0(x)$.
+as initial condition. Due to the periodic boundary conditions and the transport velocity ``a=1``, the solution at time ``t=1`` is identical to the initial distribution, i.e. ``u(1,x) = u_0(x)``.
 
 ```@setup LinearAdvection
 import Pkg; Pkg.add("PositiveIntegrators"); Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots")

From afc405f78546048780142a9704ea5ca09be38d13 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Fri, 12 Jul 2024 17:34:31 +0200
Subject: [PATCH 06/17] extended tutorial

---
 docs/src/linear_advection.md | 24 +++++++++++++++++++++---
 1 file changed, 21 insertions(+), 3 deletions(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 9983b48b..ab308a66 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -52,12 +52,12 @@ u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
 as initial condition. Due to the periodic boundary conditions and the transport velocity ``a=1``, the solution at time ``t=1`` is identical to the initial distribution, i.e. ``u(1,x) = u_0(x)``.
 
 ```@setup LinearAdvection
-import Pkg; Pkg.add("PositiveIntegrators"); Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots")
+import Pkg; Pkg.add("PositiveIntegrators"); Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots"); Pkg.add("BenchmarkTools")
 ```
 ```@example LinearAdvection
-N = 1000 # number of subintervals
+N = 100 # number of subintervals
 dx = 1/N; # mesh width
-x = LinRange(dx, 1.0, N+1) # discretization points x_1,...,x_N = x_0
+x = LinRange(dx, 1.0, N) # discretization points x_1,...,x_N = x_0
 u0 = 0.0 .+ (0.4 .≤ x .≤ 0.6) .* 1.0 # initial solution
 tspan = (0.0, 1.0) # time domain
 
@@ -98,3 +98,21 @@ plot(x,u0)
 plot!(x, last(sol.u))
 ```
 
+```@example LinearAdvection
+p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),
+                                  N - 1 => ones(eltype(u0), 1))
+prob_sparse = ConservativePDSProblem(lin_adv_P!, u0, tspan; p_prototype=p_prototype)
+
+sol_sparse = solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
+```
+
+```@example LinearAdvection
+plot(x,u0)
+plot!(x, last(sol_sparse.u))
+```
+
+```@example LinearAdvection
+using BenchmarkTools
+@benchmark solve(prob, MPRK43I(1.0, 0.5); save_everystep = false)
+@benchmark solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
+```
\ No newline at end of file

From a5f2b54ee815487bbc99ce5ab13ce781ff823d0f Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 08:15:30 +0200
Subject: [PATCH 07/17] Update docs/make.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/make.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/make.jl b/docs/make.jl
index 467a19d8..e4307ce1 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -73,7 +73,9 @@ makedocs(modules = [PositiveIntegrators],
          # Explicitly specify documentation structure
          pages = [
              "Home" => "index.md",
-             "Tutorials" => "linear_advection.md",
+             "Tutorials" => [
+                 "Linear Advection" => "linear_advection.md",
+             ],
              "API reference" => "api_reference.md",
              "Contributing" => "contributing.md",
              "Code of conduct" => "code_of_conduct.md",

From c887c57c9a3869c4d960ca67245a89b6712ad5c0 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 08:15:39 +0200
Subject: [PATCH 08/17] Update docs/src/linear_advection.md

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/src/linear_advection.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index ab308a66..290f4ae0 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -100,7 +100,7 @@ plot!(x, last(sol.u))
 
 ```@example LinearAdvection
 p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),
-                                  N - 1 => ones(eltype(u0), 1))
+                      N - 1 => ones(eltype(u0), 1))
 prob_sparse = ConservativePDSProblem(lin_adv_P!, u0, tspan; p_prototype=p_prototype)
 
 sol_sparse = solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)

From bd0bd0f3a5c8b26d09720888c1f95e81282a543c Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 08:17:04 +0200
Subject: [PATCH 09/17] Update docs/src/linear_advection.md

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/src/linear_advection.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 290f4ae0..90043a30 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -58,7 +58,7 @@ import Pkg; Pkg.add("PositiveIntegrators"); Pkg.add("OrdinaryDiffEq");  Pkg.add(
 N = 100 # number of subintervals
 dx = 1/N; # mesh width
 x = LinRange(dx, 1.0, N) # discretization points x_1,...,x_N = x_0
-u0 = 0.0 .+ (0.4 .≤ x .≤ 0.6) .* 1.0 # initial solution
+u0 = @. 0.0 + (0.4 ≤ x ≤ 0.6) * 1.0 # initial solution
 tspan = (0.0, 1.0) # time domain
 
 nothing #hide

From bf8b3bee81c7996267625c79c0f3f068dc3b8261 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 08:17:28 +0200
Subject: [PATCH 10/17] Update docs/src/linear_advection.md

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/src/linear_advection.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 90043a30..5b7796bb 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -99,6 +99,7 @@ plot!(x, last(sol.u))
 ```
 
 ```@example LinearAdvection
+using SparseArrays
 p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),
                       N - 1 => ones(eltype(u0), 1))
 prob_sparse = ConservativePDSProblem(lin_adv_P!, u0, tspan; p_prototype=p_prototype)

From a25f9b448ad4e2a32b922593edfdacde8043c27f Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Mon, 15 Jul 2024 08:43:04 +0200
Subject: [PATCH 11/17] Updated Project.toml and improved tutorial for linear
 advection

---
 docs/Project.toml            | 7 ++++++-
 docs/src/linear_advection.md | 8 +++++---
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index e7b9a6db..9dce2f22 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,9 +1,14 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PositiveIntegrators = "d1b20bf0-b083-4985-a874-dc5121669aa5"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
+BenchmarkTools = "1"
 Documenter = "1"
-OrdinaryDiffEq = "6"
+OrdinaryDiffEq = "6.59"
 Plots = "1"
+SparseArrays = "1.7"
diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 5b7796bb..064edeba 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -51,9 +51,6 @@ u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
 
 as initial condition. Due to the periodic boundary conditions and the transport velocity ``a=1``, the solution at time ``t=1`` is identical to the initial distribution, i.e. ``u(1,x) = u_0(x)``.
 
-```@setup LinearAdvection
-import Pkg; Pkg.add("PositiveIntegrators"); Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots"); Pkg.add("BenchmarkTools")
-```
 ```@example LinearAdvection
 N = 100 # number of subintervals
 dx = 1/N; # mesh width
@@ -105,6 +102,8 @@ p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),
 prob_sparse = ConservativePDSProblem(lin_adv_P!, u0, tspan; p_prototype=p_prototype)
 
 sol_sparse = solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
+
+nothing #hide
 ```
 
 ```@example LinearAdvection
@@ -115,5 +114,8 @@ plot!(x, last(sol_sparse.u))
 ```@example LinearAdvection
 using BenchmarkTools
 @benchmark solve(prob, MPRK43I(1.0, 0.5); save_everystep = false)
+```
+
+```@example LinearAdvection
 @benchmark solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
 ```
\ No newline at end of file

From 18a4a0f7a6277969497379aec385036c39b985a4 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Mon, 15 Jul 2024 09:10:51 +0200
Subject: [PATCH 12/17] N=1000 in linear advection tutorial

---
 docs/src/linear_advection.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 064edeba..ade24f15 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -52,8 +52,8 @@ u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
 as initial condition. Due to the periodic boundary conditions and the transport velocity ``a=1``, the solution at time ``t=1`` is identical to the initial distribution, i.e. ``u(1,x) = u_0(x)``.
 
 ```@example LinearAdvection
-N = 100 # number of subintervals
-dx = 1/N; # mesh width
+N = 1000 # number of subintervals
+dx = 1/N # mesh width
 x = LinRange(dx, 1.0, N) # discretization points x_1,...,x_N = x_0
 u0 = @. 0.0 + (0.4 ≤ x ≤ 0.6) * 1.0 # initial solution
 tspan = (0.0, 1.0) # time domain

From e4ab8170806569d6682ddc8eb1eaab617de6947f Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 17:56:02 +0200
Subject: [PATCH 13/17] Update docs/src/linear_advection.md

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/src/linear_advection.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index ade24f15..7e58c4bd 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -28,7 +28,7 @@ This system can also be written as ``\partial_t \mathbf u(t)=\mathbf A\mathbf u(
 \mathbf A= \frac{a}{Δ x}\begin{bmatrix}-1&0&\dots&0&1\\1&-1&\ddots&&0\\0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\0&\dots&0&1&-1\end{bmatrix}.
 ```
 
-In particular the matrix ``\mathbf A`` shows, that there is a single production term and a single destruction term per equation. 
+In particular the matrix ``\mathbf A`` shows that there is a single production term and a single destruction term per equation. 
 Furthermore, the system is conservative as ``\mathbf A`` has column sum zero.
 To be precise, the production matrix ``\mathbf P = (p_{i,j})`` of this conservative PDS is given by
 

From f3a90c9a1dbb7d7a3b8eb17d524e052528ac1131 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 17:57:10 +0200
Subject: [PATCH 14/17] Update docs/src/linear_advection.md

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/src/linear_advection.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 7e58c4bd..77ae1bb4 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -91,8 +91,8 @@ nothing #hide
 ```@example LinearAdvection
 using Plots
 
-plot(x,u0)
-plot!(x, last(sol.u))
+plot(x, u0; label = "u0", xguide = "x", yguide = "u")
+plot!(x, last(sol.u); label = "u")
 ```
 
 ```@example LinearAdvection

From 8cdb60013df4c458b8c686935113e6508662763a Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 18:02:59 +0200
Subject: [PATCH 15/17] Update docs/src/linear_advection.md

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/src/linear_advection.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 77ae1bb4..04da087b 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -95,6 +95,10 @@ plot(x, u0; label = "u0", xguide = "x", yguide = "u")
 plot!(x, last(sol.u); label = "u")
 ```
 
+### Using sparse matrices
+
+TODO: Some text
+
 ```@example LinearAdvection
 using SparseArrays
 p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),

From 4d845d71a1e216eb656eb610d87bd2b2453dfb3f Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <SKopecz@users.noreply.github.com>
Date: Mon, 15 Jul 2024 18:03:13 +0200
Subject: [PATCH 16/17] Update docs/src/linear_advection.md

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 docs/src/linear_advection.md | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 04da087b..70f2d0e4 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -115,6 +115,11 @@ plot(x,u0)
 plot!(x, last(sol_sparse.u))
 ```
 
+### Performance comparison
+
+Finally, we use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl)
+to compare the performance of the different implementations.
+
 ```@example LinearAdvection
 using BenchmarkTools
 @benchmark solve(prob, MPRK43I(1.0, 0.5); save_everystep = false)

From c37dbba9073efb335cf0b6795fad2ec1840a4bd5 Mon Sep 17 00:00:00 2001
From: Stefan Kopecz <kopecz@mathematik.uni-kassel.de>
Date: Mon, 15 Jul 2024 18:23:13 +0200
Subject: [PATCH 17/17] Added plot labels

---
 docs/src/linear_advection.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 70f2d0e4..6694d441 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -111,8 +111,8 @@ nothing #hide
 ```
 
 ```@example LinearAdvection
-plot(x,u0)
-plot!(x, last(sol_sparse.u))
+plot(x,u0; label = "u0", xguide = "x", yguide = "u")
+plot!(x, last(sol_sparse.u); label = "u")
 ```
 
 ### Performance comparison