diff --git a/.github/workflows/DocPreviewCleanup.yml b/.github/workflows/DocPreviewCleanup.yml
index 7bd74fac3..02eddbcf6 100644
--- a/.github/workflows/DocPreviewCleanup.yml
+++ b/.github/workflows/DocPreviewCleanup.yml
@@ -8,7 +8,7 @@ jobs:
   doc-preview-cleanup:
     # Do not run on forks to avoid authorization errors
     # Source: https://github.community/t/have-github-action-only-run-on-master-repo-and-not-on-forks/140840/18
-    if: github.repository_owner == 'SKopecz'
+    if: github.repository_owner == 'NumericalMathematics'
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
diff --git a/CITATION.bib b/CITATION.bib
index 19c2a16a6..bd7352c12 100644
--- a/CITATION.bib
+++ b/CITATION.bib
@@ -4,5 +4,5 @@ @misc{PositiveIntegrators.jl
   author={Kopecz, Stefan and Lampert, Joshua and Ranocha, Hendrik},
   year={2023},
   doi={10.5281/zenodo.10868393},
-  url={https://github.com/SKopecz/PositiveIntegrators.jl}
+  url={https://github.com/NumericalMathematics/PositiveIntegrators.jl}
 }
diff --git a/Project.toml b/Project.toml
index 1c90025a4..77fa65056 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PositiveIntegrators"
 uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
 authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
-version = "0.2.11-DEV"
+version = "0.2.11"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
diff --git a/README.md b/README.md
index cad0150d9..f00f9d163 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,10 @@
 # PositiveIntegrators.jl
 
-[![Docs-stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://SKopecz.github.io/PositiveIntegrators.jl/stable)
-[![Docs-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://SKopecz.github.io/PositiveIntegrators.jl/dev)
-[![Build Status](https://github.com/SKopecz/PositiveIntegrators.jl/workflows/CI/badge.svg)](https://github.com/SKopecz/PositiveIntegrators.jl/actions?query=workflow%3ACI)
-[![Codecov](https://codecov.io/gh/SKopecz/PositiveIntegrators.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/SKopecz/PositiveIntegrators.jl)
-[![Coveralls](https://coveralls.io/repos/github/SKopecz/PositiveIntegrators.jl/badge.svg?branch=main)](https://coveralls.io/github/SKopecz/PositiveIntegrators.jl?branch=main)
+[![Docs-stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://NumericalMathematics.github.io/PositiveIntegrators.jl/stable)
+[![Docs-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://NumericalMathematics.github.io/PositiveIntegrators.jl/dev)
+[![Build Status](https://github.com/NumericalMathematics/PositiveIntegrators.jl/workflows/CI/badge.svg)](https://github.com/NumericalMathematics/PositiveIntegrators.jl/actions?query=workflow%3ACI)
+[![Codecov](https://codecov.io/gh/NumericalMathematics/PositiveIntegrators.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/NumericalMathematics/PositiveIntegrators.jl)
+[![Coveralls](https://coveralls.io/repos/github/NumericalMathematics/PositiveIntegrators.jl/badge.svg?branch=main)](https://coveralls.io/github/NumericalMathematics/PositiveIntegrators.jl?branch=main)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10868393.svg)](https://doi.org/10.5281/zenodo.10868393)
@@ -20,12 +20,12 @@ by
   modified Patankar-Runge-Kutta (MPRK) schemes
 
 Please find more information online in the
-[documentation](https://skopecz.github.io/PositiveIntegrators.jl/stable).
+[documentation](https://NumericalMathematics.github.io/PositiveIntegrators.jl/stable).
 
 
 ## Installation
 
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl)
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl)
 is a registered Julia package. Thus, you can install it from the Julia REPL via
 ```julia
 julia> using Pkg; Pkg.add("PositiveIntegrators")
@@ -40,7 +40,7 @@ julia> using Pkg; Pkg.update("PositiveIntegrators")
 ## Referencing
 
 If you use
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl)
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl)
 for your research, please cite it using the bibtex entry
 ```bibtex
 @misc{PositiveIntegrators.jl,
@@ -49,7 +49,7 @@ for your research, please cite it using the bibtex entry
   author={Kopecz, Stefan and Lampert, Joshua and Ranocha, Hendrik},
   year={2023},
   doi={10.5281/zenodo.10868393},
-  url={https://github.com/SKopecz/PositiveIntegrators.jl}
+  url={https://github.com/NumericalMathematics/PositiveIntegrators.jl}
 }
 ```
 
@@ -57,8 +57,8 @@ for your research, please cite it using the bibtex entry
 ## License and contributing
 
 This project is licensed under the MIT license
-(see [License](https://github.com/SKopecz/PositiveIntegrators.jl/blob/main/LICENSE)).
+(see [License](https://github.com/NumericalMathematics/PositiveIntegrators.jl/blob/main/LICENSE)).
 Since it is an open-source project, we are very happy to accept contributions
 from the community. Please refer to the section
-[Contributing](https://github.com/SKopecz/PositiveIntegrators.jl/blob/main/CONTRIBUTING.md)
+[Contributing](https://github.com/NumericalMathematics/PositiveIntegrators.jl/blob/main/CONTRIBUTING.md)
 for more details.
diff --git a/SECURITY.md b/SECURITY.md
index 8336024eb..48b5e886d 100644
--- a/SECURITY.md
+++ b/SECURITY.md
@@ -14,7 +14,7 @@ used in the Julia ecosystem is supported with security updates.
 ## Reporting a Vulnerability
 
 To report a security issue, please use the GitHub Security Advisory
-["Report a Vulnerability"](https://github.com/SKopecz/PositiveIntegrators.jl/security/advisories/new)
+["Report a Vulnerability"](https://github.com/NumericalMathematics/PositiveIntegrators.jl/security/advisories/new)
 tab.
 
 We will send a response indicating the next steps in handling your report.
diff --git a/docs/make.jl b/docs/make.jl
index d54c7bbc4..9a4786946 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -18,11 +18,12 @@ DocMeta.setdocmeta!(PositiveIntegrators,
 # as necessary
 open(joinpath(@__DIR__, "src", "license.md"), "w") do io
     # Point to source license file
-    println(io, """
-    ```@meta
-    EditURL = "https://github.com/SKopecz/PositiveIntegrators.jl/blob/main/LICENSE"
-    ```
-    """)
+    println(io,
+            """
+```@meta
+EditURL = "https://github.com/NumericalMathematics/PositiveIntegrators.jl/blob/main/LICENSE"
+```
+""")
     # Write the modified contents
     println(io, "# License")
     println(io, "")
@@ -37,7 +38,7 @@ open(joinpath(@__DIR__, "src", "code_of_conduct.md"), "w") do io
     println(io,
             """
             ```@meta
-            EditURL = "https://github.com/SKopecz/PositiveIntegrators.jl/blob/main/CODE_OF_CONDUCT.md"
+            EditURL = "https://github.com/NumericalMathematics/PositiveIntegrators.jl/blob/main/CODE_OF_CONDUCT.md"
             ```
             """)
     # Write the modified contents
@@ -51,11 +52,12 @@ end
 
 open(joinpath(@__DIR__, "src", "contributing.md"), "w") do io
     # Point to source license file
-    println(io, """
-    ```@meta
-    EditURL = "https://github.com/SKopecz/PositiveIntegrators.jl/blob/main/CONTRIBUTING.md"
-    ```
-    """)
+    println(io,
+            """
+```@meta
+EditURL = "https://github.com/NumericalMathematics/PositiveIntegrators.jl/blob/main/CONTRIBUTING.md"
+```
+""")
     # Write the modified contents
     println(io, "# Contributing")
     println(io, "")
@@ -69,7 +71,7 @@ end
 makedocs(modules = [PositiveIntegrators],
          sitename = "PositiveIntegrators.jl",
          format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
-                                  canonical = "https://SKopecz.github.io/PositiveIntegrators.jl/stable"),
+                                  canonical = "https://NumericalMathematics.github.io/PositiveIntegrators.jl/stable"),
          # Explicitly specify documentation structure
          pages = [
              "Home" => "index.md",
@@ -95,6 +97,6 @@ makedocs(modules = [PositiveIntegrators],
              "License" => "license.md"
          ])
 
-deploydocs(repo = "github.com/SKopecz/PositiveIntegrators.jl",
+deploydocs(repo = "github.com/NumericalMathematics/PositiveIntegrators.jl",
            devbranch = "main",
            push_preview = true)
diff --git a/docs/src/faq.md b/docs/src/faq.md
index a112682bb..36dcf683d 100644
--- a/docs/src/faq.md
+++ b/docs/src/faq.md
@@ -3,7 +3,7 @@
 ## Sparse matrices
 
 You can use sparse matrices for the linear systems arising in
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl),
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl),
 as described, e.g., in the [tutorial on linear advection](@ref tutorial-linear-advection).
 However, you need to make sure that you do not change the sparsity pattern
 of the production term matrix since we assume that the structural nonzeros
diff --git a/docs/src/heat_equation_dirichlet.md b/docs/src/heat_equation_dirichlet.md
index 8b1d88ffb..9569ac840 100644
--- a/docs/src/heat_equation_dirichlet.md
+++ b/docs/src/heat_equation_dirichlet.md
@@ -58,7 +58,7 @@ In addition, all production and destruction terms not listed are zero.
 
 Now we are ready to define a [`PDSProblem`](@ref) and to solve this
 problem with a method of
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) or
 [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 In the following we use ``N = 100`` nodes and the time domain ``t \in [0,1]``.
 Moreover, we choose the initial condition
diff --git a/docs/src/heat_equation_neumann.md b/docs/src/heat_equation_neumann.md
index a1b6105ce..446f10cfe 100644
--- a/docs/src/heat_equation_neumann.md
+++ b/docs/src/heat_equation_neumann.md
@@ -50,7 +50,7 @@ and destruction terms ``d_{i,j} = p_{j,i}``. In addition, all production and des
 
 Now we are ready to define a [`ConservativePDSProblem`](@ref) and to solve this
 problem with a method of
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) or
 [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 In the following we use ``N = 100`` nodes and the time domain ``t \in [0,1]``.
 Moreover, we choose the initial condition
diff --git a/docs/src/index.md b/docs/src/index.md
index 2021fdf51..4fcc56a4c 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,14 +1,14 @@
 # PositiveIntegrators.jl
 
 The [Julia](https://julialang.org) library
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl)
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl)
 provides several time integration methods developed to preserve the positivity
 of numerical solutions.
 
 
 ## Installation
 
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl)
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl)
 is a registered Julia package. Thus, you can install it from the Julia REPL via
 ```julia
 julia> using Pkg; Pkg.add("PositiveIntegrators")
@@ -88,7 +88,7 @@ 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 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.
+Now that the problem has been created, we can solve it with any method of [PositiveIntegrators.jl](https://github.com/NumericalMathematics/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
 sol = solve(prob, MPRK22(1.0))
@@ -132,7 +132,7 @@ The corresponding production matrix ``\mathbf P`` is
 \mathbf P(S,I,R) = \begin{pmatrix}0 & 0 & 0\\ \frac{β S I}{N} & 0 & 0\\ 0 & γ I & 0\end{pmatrix}.
 ```
 
-The following example shows how to implement the above SIR model with ``\beta=0.4, \gamma=0.04``, initial conditions ``S(0)=997, I(0)=3, R(0)=0`` and time domain ``(0, 100)`` using `ConservativePDSProblem` from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+The following example shows how to implement the above SIR model with ``\beta=0.4, \gamma=0.04``, initial conditions ``S(0)=997, I(0)=3, R(0)=0`` and time domain ``(0, 100)`` using `ConservativePDSProblem` from [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl).
 
 ```@example SIR
 using PositiveIntegrators
@@ -158,7 +158,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 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/NumericalMathematics/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
@@ -186,7 +186,7 @@ for your research, please cite it using the bibtex entry
   author={Kopecz, Stefan and Lampert, Joshua and Ranocha, Hendrik},
   year={2023},
   doi={10.5281/zenodo.10868393},
-  url={https://github.com/SKopecz/PositiveIntegrators.jl}
+  url={https://github.com/NumericalMathematics/PositiveIntegrators.jl}
 }
 ```
 
diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md
index 457bc7618..ab78f435d 100644
--- a/docs/src/linear_advection.md
+++ b/docs/src/linear_advection.md
@@ -42,7 +42,7 @@ In addition, all production and destruction terms not listed have the value zero
 
 ## Solution of the production-destruction system
 
-Now we are ready to define a [`ConservativePDSProblem`](@ref) 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`](@ref) and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/NumericalMathematics/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}
diff --git a/docs/src/npzd_model.md b/docs/src/npzd_model.md
index 9230ff9c6..16184d279 100644
--- a/docs/src/npzd_model.md
+++ b/docs/src/npzd_model.md
@@ -28,7 +28,7 @@ whereby production terms not listed have the value zero. Since the PDS is conser
 
 ## Solution of the production-destruction system
 
-Now we are ready to define a [`ConservativePDSProblem`](@ref) 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/).
+Now we are ready to define a [`ConservativePDSProblem`](@ref) and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 
 As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are
 1. an out-of-place implementation with standard (dynamic) matrices and vectors,
@@ -76,7 +76,7 @@ using Plots
 
 plot(sol_oop; label = ["N" "P" "Z" "D"], xguide = "t")
 ```
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) provides the function [`isnonnegative`](@ref) (and also [`isnegative`](@ref)) to check if the solution is actually nonnegative, as expected from an MPRK scheme.
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) provides the function [`isnonnegative`](@ref) (and also [`isnegative`](@ref)) to check if the solution is actually nonnegative, as expected from an MPRK scheme.
 ```@example NPZD
 isnonnegative(sol_oop)
 ```
@@ -173,7 +173,7 @@ This solution is also nonnegative.
 isnonnegative(sol_static)
 ```
 
-The above implementation of the NPZD model using `StaticArrays` can also be found in the [Example Problems](https://skopecz.github.io/PositiveIntegrators.jl/dev/api_reference/#Example-problems) as [`prob_pds_npzd`](@ref).
+The above implementation of the NPZD model using `StaticArrays` can also be found in the [Example Problems](https://NumericalMathematics.github.io/PositiveIntegrators.jl/dev/api_reference/#Example-problems) as [`prob_pds_npzd`](@ref).
 
 ### Performance comparison
 
diff --git a/docs/src/npzd_model_benchmark.md b/docs/src/npzd_model_benchmark.md
index b3e09f232..eda30bfcb 100644
--- a/docs/src/npzd_model_benchmark.md
+++ b/docs/src/npzd_model_benchmark.md
@@ -1,6 +1,6 @@
 # [Benchmark: Solution of an NPZD model](@id benchmark-npzd)
 
-We use the NPZD model [`prob_pds_npzd`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+We use the NPZD model [`prob_pds_npzd`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl).
 
 ```@example NPZD
 using OrdinaryDiffEqLowOrderRK, OrdinaryDiffEqSDIRK, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqTsit5, OrdinaryDiffEqVerner
@@ -19,14 +19,14 @@ using Plots
 function npzd_plot(sol, sol_ref = nothing, title = "")
      colors = palette(:default)[1:4]'
      if !isnothing(sol_ref)
-          p = plot(sol_ref, linestyle = :dash, label = "", color = colors, 
+          p = plot(sol_ref, linestyle = :dash, label = "", color = colors,
                    linewidth = 2)
           plot!(p, sol; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
                 color = colors, title, label = ["N" "P" "Z" "D"], legend = :right,
                 linewidth = 2);
      else
           p = plot(sol; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
-                   color = colors, title, label = ["N" "P" "Z" "D"], legend = :right, 
+                   color = colors, title, label = ["N" "P" "Z" "D"], legend = :right,
                    linewidths = 2);
      end
      return p
@@ -35,8 +35,8 @@ nothing  # hide
 ```
 
 Standard methods have difficulties to solve the NPZD problem accurately for loose tolerances or large time step sizes.
-This is because the first variable, ``N``, has only a tiny margin for negative values. 
-In most cases, negative values of ``N``will directly decrease ``N``further, resulting in completely inaccurate solutions.  
+This is because the first variable, ``N``, has only a tiny margin for negative values.
+In most cases, negative values of ``N``will directly decrease ``N``further, resulting in completely inaccurate solutions.
 
 ```@example NPZD
 # compute reference solution for plotting
@@ -50,14 +50,14 @@ sol_MPRK = solve(prob, MPRK22(1.0); abstol, reltol);
 
 # plot solutions
 p1 = npzd_plot(sol_Ros23, ref_sol, "Rosenbrock23"); # helper function defined above
-p2 = npzd_plot(sol_MPRK, ref_sol, "MPRK22(1.0)"); 
+p2 = npzd_plot(sol_MPRK, ref_sol, "MPRK22(1.0)");
 plot(p1, p2)
 ```
 
-Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions. 
+Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions.
 
 ```@example NPZD
-sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol, 
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol,
                   isoutofdomain = isnegative); #reject negative solutions
 
 npzd_plot(sol_Ros23, ref_sol) #auxiliary function defined above
@@ -101,9 +101,9 @@ We start with a comparison of different adaptive MPRK schemes.
 
 ```@example NPZD
 # choose methods to compare
-algs = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0); 
+algs = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0);
         MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0)]
-labels = ["MPRK22(0.5)"; "MPPRK22(2/3)"; "MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)"; 
+labels = ["MPRK22(0.5)"; "MPPRK22(2/3)"; "MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)";
           "MPRK43I(1.0, 0.5)"; "MPRK43I(0.5, 0.75)"; "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)"]
 
 # compute work-precision data
@@ -111,7 +111,7 @@ wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
                                compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "NPZD benchmark", legend = :topright,     
+plot(wp, labels; title = "NPZD benchmark", legend = :topright,
      color = permutedims([repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...]),
      xlims = (10^-7, 2*10^-1), xticks = 10.0 .^ (-8:1:0),
      ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
@@ -128,7 +128,7 @@ p2 = npzd_plot(sol_MPRK43, ref_sol, "MPRK43I(1.0, 0.5)");
 plot(p1, p2)
 ```
 
-Next we compare `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)` to explicit and implicit methods of second and third order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). 
+Next we compare `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)` to explicit and implicit methods of second and third order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 To guarantee nonnegative solutions, we select the solver option `isoutofdomain = isnegative`.
 
 ```@example NPZD
@@ -160,14 +160,14 @@ We see that for the NPZD problem the use of adaptive MPRK schemes is only benefi
 Now we compare `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)` to [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Again, to guarantee positive solutions we select the solver option `isoutofdomain = isnegative`.
 
 ```@example NPZD
-algs3 = [Tsit5(); BS3(); Vern6(); Vern7(); Vern8(); TRBDF2(); Rosenbrock23(); 
+algs3 = [Tsit5(); BS3(); Vern6(); Vern7(); Vern8(); TRBDF2(); Rosenbrock23();
          Rodas5P(); Rodas4P()]
 labels3 = ["Tsit5"; "BS3"; "Vern6"; "Vern7"; "Vern8"; "TRBDF2"; "Rosenbrock23";
           "Rodas5P"; "Rodas4P"]
 
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
-                               compute_error) 
+                               compute_error)
 # add work-precision data with isoutofdomain = isnegative
 work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
                                compute_error, isoutofdomain = isnegative)
@@ -183,7 +183,7 @@ We see that it is advisable to use a high-order explicit method like `Vern7()` o
 
 #### Relative maximum error over all time steps
 
-In this section we do not compare the relative maximum errors at the final time ``t = 10.0``, but the relative maximum errors over all time steps. 
+In this section we do not compare the relative maximum errors at the final time ``t = 10.0``, but the relative maximum errors over all time steps.
 
 ```@example NPZD
 # select relative maximum error over all time steps
@@ -191,8 +191,8 @@ compute_error = rel_max_error_overall
 nothing  # hide
 ```
 
-The results are very similar to those from above. 
-We therefore only show the work-precision diagrams without further comments. 
+The results are very similar to those from above.
+We therefore only show the work-precision diagrams without further comments.
 The main difference are significantly increased errors which mainly occur around time ``t = 2.0`` where there is a sharp kink in the solution.
 
 ```@example NPZD
@@ -201,7 +201,7 @@ wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
                                compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "NPZD benchmark", legend = :topright,     
+plot(wp, labels; title = "NPZD benchmark", legend = :topright,
           color = permutedims([repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...]),
           xlims = (10^-5, 10^4), xticks = 10.0 .^ (-5:1:4),
           ylims = (10^-6, 10^-1), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
@@ -226,7 +226,7 @@ plot(wp, [labels1; labels2]; title = "NPZD benchmark", legend = :topright,
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
                                compute_error)
-# add work-precision data with isoutofdomain = isnegative                             
+# add work-precision data with isoutofdomain = isnegative
 work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
                                compute_error, isoutofdomain=isnegative)
 
@@ -239,8 +239,8 @@ plot(wp, [labels1; labels3]; title = "NPZD benchmark", legend = :topright,
 
 ### Fixed time step sizes
 
-Here we use fixed time step sizes instead of adaptive time stepping. 
-Similar to the adaptive situation above, standard schemes are likely to compute negative solutions for the NPZD problem. 
+Here we use fixed time step sizes instead of adaptive time stepping.
+Similar to the adaptive situation above, standard schemes are likely to compute negative solutions for the NPZD problem.
 
 ```@example NPZD
 sol_Ros23 = solve(prob, Rosenbrock23(), dt = 1.0, adaptive = false);
@@ -252,7 +252,7 @@ plot(p1, p2)
 ```
 
 We use the functions [`work_precision_fixed`](@ref) and [`work_precision_fixed!`](@ref) to compute the data for the diagrams.
-Please note that these functions set error and computing time to `Inf`, whenever a solution contains negative elements. 
+Please note that these functions set error and computing time to `Inf`, whenever a solution contains negative elements.
 Consequently, such cases are not visible in the work-precision diagrams.
 
 Within the work-precision diagrams we use the following time step sizes.
@@ -273,7 +273,7 @@ compute_error = rel_max_error_tend
 nothing  # hide
 ```
 
-First, we compare different MPRK methods. 
+First, we compare different MPRK methods.
 For fixed time step sizes we can also consider `MPE()` and `SSPMPRK43()`.
 
 ```@example NPZD
@@ -286,13 +286,13 @@ wp = work_precision_fixed(prob, algs, labels, dts, alg_ref;
                           compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "NPZD benchmark", legend = :bottomleft,     
+plot(wp, labels; title = "NPZD benchmark", legend = :bottomleft,
      color = permutedims([5,repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...,6]),
      xlims = (10^-10, 1*10^0), xticks = 10.0 .^ (-10:1:0),
      ylims = (1*10^-6, 10^-1), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
 ```
 
-Apart from `MPE()` the schemes behave very similar and a difference in order can only be observed for the smaller step sizes. 
+Apart from `MPE()` the schemes behave very similar and a difference in order can only be observed for the smaller step sizes.
 We choose `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)` for comparisons with other second- and third-order schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 
 ```@example NPZD
@@ -304,7 +304,7 @@ wp = work_precision_fixed(prob, [algs1; algs2], [labels1; labels2], dts, alg_ref
 plot(wp, [labels1; labels2]; title = "NPZD benchmark", legend = :topright,
      color = permutedims([1, 3, repeat([4], 3)..., repeat([5],4)...,repeat([6],4)...]),
      xlims = (10^-13, 10^2), xticks = 10.0 .^ (-12:2:6),
-     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
 ```
 
 We see that the MPRK schemes are to be preferred for the rather large step sizes ``\Delta t \in\lbrace 1.0, 0.5, 0.25, 0.125\rbrace``, for which the other schemes cannot provide nonnegative solutions.
@@ -328,7 +328,7 @@ wp = work_precision_fixed(prob, [algs1; algs3], [labels1; labels3], dts, alg_ref
 plot(wp, [labels1; labels3]; title = "NPZD benchmark", legend = :topright,
      color = permutedims([1, 3, repeat([4], 3)..., repeat([5],4)...,repeat([6],4)...]),
      xlims = (10^-14, 10^0), xticks = 10.0 .^ (-14:2:10),
-     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
 ```
 
 #### Relative maximum error over all time steps
@@ -348,7 +348,7 @@ wp = work_precision_fixed(prob, algs, labels, dts, alg_ref;
                                compute_error)
 
 #plot work-precision diagram
-plot(wp, labels; title = "NPZD benchmark", legend = :bottomleft,     
+plot(wp, labels; title = "NPZD benchmark", legend = :bottomleft,
      color = permutedims([5,repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...,6]),
      xlims = (10^-4, 10^5), xticks = 10.0 .^ (-4:1:5),
      ylims = (10^-6, 10^-1), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
@@ -358,24 +358,24 @@ plot(wp, labels; title = "NPZD benchmark", legend = :bottomleft,
 wp = work_precision_fixed(prob, algs1, labels1, dts, alg_ref;
                                compute_error)
 work_precision_fixed!(wp, prob, algs2, labels2, dts, alg_ref;
-                     compute_error)                               
+                     compute_error)
 
 plot(wp, [labels1; labels2]; title = "NPZD benchmark", legend = :topright,
      color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
      xlims = (10^-4, 10^6), xticks = 10.0 .^ (-12:1:6),
-     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
 ```
 
 ```@example NPZD
 wp = work_precision_fixed(prob, algs1, labels1, dts, alg_ref;
                                compute_error)
 work_precision_fixed!(wp, prob, algs3, labels3, dts, alg_ref;
-                     compute_error)                               
+                     compute_error)
 
 plot(wp, [labels1; labels3]; title = "NPZD benchmark", legend = :bottomleft,
      color = permutedims([1, 3, repeat([4], 5)..., 5, repeat([7], 3)...]),
      xlims = (10^-12, 10^6), xticks = 10.0 .^ (-12:2:6),
-     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
 ```
 
 ## Package versions
diff --git a/docs/src/robertson.md b/docs/src/robertson.md
index 3cd1d7e76..381c5cda8 100644
--- a/docs/src/robertson.md
+++ b/docs/src/robertson.md
@@ -27,7 +27,7 @@ whereby production terms not listed have the value zero. Since the PDS is conser
 
 ## Solution of the production-destruction system
 
-Now we are ready to define a [`ConservativePDSProblem`](@ref) and to solve this problem with any method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) which is suited for stiff problems.
+Now we are ready to define a [`ConservativePDSProblem`](@ref) and to solve this problem with any method of [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) which is suited for stiff problems.
 
 Since this PDS consists of only three differential equations we provide an out-of-place implementation for the production matrix. Furthermore, we use static arrays from [StaticArrays.jl](https://juliaarrays.github.io/StaticArrays.jl/stable/) for additional efficiency. See also the tutorials on the solution of [an NPZD model](@ref tutorial-npzd) or [an stratospheric reaction problem](@ref tutorial-stratos).
 
@@ -53,7 +53,7 @@ plot(sol, tspan = (1e-6, 1e11),  xaxis = :log,
      idxs = [(0, 1), ((x, y) -> (x, 1e4 .* y), 0, 2), (0, 3)],
      label = ["u₁" "10⁴u₂" "u₃"])
 ```
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) provides the function [`isnonnegative`](@ref) (and also [`isnegative`](@ref)) to check if the solution is actually nonnegative, as expected from an MPRK scheme.
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) provides the function [`isnonnegative`](@ref) (and also [`isnegative`](@ref)) to check if the solution is actually nonnegative, as expected from an MPRK scheme.
 ```@example robertson
 isnonnegative(sol)
 ```
diff --git a/docs/src/robertson_benchmark.md b/docs/src/robertson_benchmark.md
index 639cdc8a9..8b84fc04e 100644
--- a/docs/src/robertson_benchmark.md
+++ b/docs/src/robertson_benchmark.md
@@ -1,6 +1,6 @@
 # [Benchmark: Solution of the Robertson problem](@id benchmark-robertson)
 
-Here we use the stiff Robertson problem [`prob_pds_robertson`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+Here we use the stiff Robertson problem [`prob_pds_robertson`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl).
 
 ```@example ROBER
 using OrdinaryDiffEqFIRK, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqSDIRK
@@ -39,7 +39,7 @@ end
 nothing # hide
 ```
 
-For this stiff problem the computation of negative approximations may lead to inaccurate solutions. 
+For this stiff problem the computation of negative approximations may lead to inaccurate solutions.
 This typically occurs when adaptive time stepping uses loose tolerances.
 
 ```@example ROBER
@@ -58,10 +58,10 @@ p2 = robertson_plot(sol_MPRK, ref_sol, "MPRK22(1.0)");
 plot(p1, p2)
 ```
 
-Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions. 
+Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions.
 
 ```@example ROBER
-sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol, 
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol,
                   isoutofdomain = isnegative) #reject negative solutions
 
 robertson_plot(sol_Ros23, ref_sol, "Rosenbrock23")
@@ -69,7 +69,7 @@ robertson_plot(sol_Ros23, ref_sol, "Rosenbrock23")
 
 ## Work-Precision diagrams
 
-In the following we show several work-precision diagrams, which compare different methods with respect to computing time and the respective error. 
+In the following we show several work-precision diagrams, which compare different methods with respect to computing time and the respective error.
 We focus solely on adaptive methods, since the time interval ``(0, 10^{11})`` is too large to generate accurate solutions with fixed step sizes.
 
 Since the Robertson problem is stiff, we need to use a suited implicit scheme to compute a reference solution, see the [solver guide](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/#Stiff-Problems).
@@ -114,7 +114,7 @@ wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
                             adaptive_ref = true, compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "Robertson benchmark", legend = :topright,     
+plot(wp, labels; title = "Robertson benchmark", legend = :topright,
      color = permutedims([repeat([1], 3)..., repeat([3], 2)..., repeat([4], 2)...]),
      xlims = (10^-10, 10^0), xticks = 10.0 .^ (-10:1:0),
      ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
@@ -159,7 +159,7 @@ labels2 = ["TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbro
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error)
-# add work-precision data with isoutofdomain = isnegative                               
+# add work-precision data with isoutofdomain = isnegative
 work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error, isoutofdomain=isnegative)
 
@@ -181,7 +181,7 @@ labels3 = ["Rodas5P"; "Rodas4P"; "RadauIIA5"]
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error)
-# add work-precision data with isoutofdomain = isnegative                                 
+# add work-precision data with isoutofdomain = isnegative
 work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error, isoutofdomain=isnegative)
 
@@ -196,7 +196,7 @@ Again, we see that the MPRK schemes are in general only beneficial if low accura
 
 ### Relative maximum error over all time steps
 
-In this section we do not compare the relative maximum errors at the final time ``t = 10^{11}``, but the relative maximum errors over all time steps. 
+In this section we do not compare the relative maximum errors at the final time ``t = 10^{11}``, but the relative maximum errors over all time steps.
 
 ```@example ROBER
 # select relative maximum error at the end of the problem's time span.
@@ -212,7 +212,7 @@ wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
                             adaptive_ref = true, compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "Robertson benchmark", legend = :top,     
+plot(wp, labels; title = "Robertson benchmark", legend = :top,
      color = permutedims([repeat([1], 3)..., repeat([3], 2)..., repeat([4], 2)...]),
      xlims = (10^-4, 5*10^1), xticks = 10.0 .^ (-5:1:2),
      ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
@@ -230,7 +230,7 @@ labels1 = ["MPRK22(1.0)"; "MPRK43I(0.5,0.75)"]
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error)
-# add work-precision data with isoutofdomain = isnegative                               
+# add work-precision data with isoutofdomain = isnegative
 work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error, isoutofdomain=isnegative)
 
@@ -249,7 +249,7 @@ Finally, we compare `MPRK43I(0.5, 0.75)` and `MPRK22(1.0)` to [recommended solve
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error)
-# add work-precision data with isoutofdomain = isnegative                                 
+# add work-precision data with isoutofdomain = isnegative
 work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
                                adaptive_ref = true, compute_error, isoutofdomain=isnegative)
 
diff --git a/docs/src/scalar_pds.md b/docs/src/scalar_pds.md
index 6ac5ed269..0cc6d2510 100644
--- a/docs/src/scalar_pds.md
+++ b/docs/src/scalar_pds.md
@@ -1,7 +1,7 @@
 # [Tutorial: Solving a scalar production-destruction equation with MPRK schemes](@id tutorial-scalar-pds)
 
 Originally, modified Patankar-Runge-Kutta (MPRK) schemes were designed to solve positive and conservative systems of ordinary differential equations.
-The conservation property requires that the system consists of at least two scalar differential equations. 
+The conservation property requires that the system consists of at least two scalar differential equations.
 Nevertheless, we can also apply the idea of the Patankar trick to a scalar production-destruction system (PDS)
 
 ```math
@@ -9,15 +9,15 @@ u'(t)=p(u(t))-d(u(t)),\quad u(0)=u_0>0
 ```
 
 with nonnegative functions ``p`` and ``d``.
-Since conservation is not an issue here, we can apply the Patankar trick to the destruction term ``d`` to ensure positivity and leave the production term ``p`` unweighted. 
+Since conservation is not an issue here, we can apply the Patankar trick to the destruction term ``d`` to ensure positivity and leave the production term ``p`` unweighted.
 A first-order scheme of this type, based on the forward Euler method, reads
 
 ```math
 u^{n+1}= u^n + Δ t p(u^n) - Δ t d(u^n)\frac{u^{n+1}}{u^n}
 ```
-and this idea can easily be generalized to higher-order explicit Runge-Kutta schemes. 
+and this idea can easily be generalized to higher-order explicit Runge-Kutta schemes.
 
-By closer inspection we realize that this is exactly the approach the MPRK schemes of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) use to solve non-conservative PDS for which the production matrix is diagonal. 
+By closer inspection we realize that this is exactly the approach the MPRK schemes of [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) use to solve non-conservative PDS for which the production matrix is diagonal.
 Hence, we can use the existing schemes to solve a scalar PDS by regarding the production term as a ``1×1``-matrix and the destruction term as a ``1``-vector.
 
 # [Example 1](@id scalar-example-1)
@@ -50,7 +50,7 @@ tspan = (0.0, 10.0)
 # Attention: Input u is a 1-vector
 prod(u, p, t) = @SMatrix [u[1]^2] # create static 1x1-matrix
 dest(u, p, t) = @SVector [u[1]] # create static 1-vector
-prob = PDSProblem(prod, dest, u0, tspan) 
+prob = PDSProblem(prod, dest, u0, tspan)
 
 sol = solve(prob, MPRK22(1.0))
 
@@ -63,27 +63,27 @@ plot!(sol, label="u")
 
 # [Example 2] (@id scalar-example-2)
 
-Next, we want to compute positive solutions of a more challenging scalar PDS. 
+Next, we want to compute positive solutions of a more challenging scalar PDS.
 In [Example 1](@ref scalar-example-1), we could have also used standard schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and use the solver option `isoutofdomain` to ensure positivity.
 But this is not always the case as the following example will show.
 
-We want to compute the nonnegative solution of 
+We want to compute the nonnegative solution of
 
 ```math
-u'(t) = -\sqrt{\lvert u(t)\rvert },\quad u(0)=1. 
+u'(t) = -\sqrt{\lvert u(t)\rvert },\quad u(0)=1.
 ```
 
 for ``t≥ 0``.
 Please note that this initial value problem has infinitely many solutions
 
-```math 
+```math
 u(t) = \begin{cases} \frac{1}{4}(t-2)^2, & 0≤ t< 2,\\ 0, & 2≤ t < t^*,\\ -\frac{1}{4}(t-2)^2, & t^*≤  t, \end{cases}
 ```
 
 where ``t^*≥ 2`` is arbitrary.
 But among these, the only nonnegative solution is
 
-```math 
+```math
 u(t) = \begin{cases} \frac{1}{4}(t-2)^2, & 0≤ t< 2,\\ 0, & 2≤ t. \end{cases}
 ```
 
@@ -104,11 +104,11 @@ prob = ODEProblem(f, u0, tspan)
 sol = solve(prob, Rosenbrock23(); isoutofdomain = (u, p, t) -> any(<(0), u))
 ```
 
-We see that `isoutofdomain` cannot be used to ensure nonnegative solutions in this case, as the computation stops at about ``t≈ 2`` before the desired final time is reached. 
+We see that `isoutofdomain` cannot be used to ensure nonnegative solutions in this case, as the computation stops at about ``t≈ 2`` before the desired final time is reached.
 For at least first- and second-order explicit Runge-Kutta schemes, this can also be shown analytically. A brief computation reveals that to ensure nonnegative solutions, the time step size must tend to zero if the numerical solution tends to zero.
 
-Next, we want to use an MPRK scheme. 
-We can choose ``p(u)=0`` as the production term and ``d(u)=\sqrt{\lvert u\rvert }`` as the destruction term. 
+Next, we want to use an MPRK scheme.
+We can choose ``p(u)=0`` as the production term and ``d(u)=\sqrt{\lvert u\rvert }`` as the destruction term.
 Furthermore, we create the [`PDSProblem`](@ref) in the same way as in [Example 1](@ref scalar-example-1).
 
 ```@example
@@ -130,7 +130,7 @@ plot(tt, f.(tt), label="exact")
 plot!(sol, label="u")
 ```
 
-We can see that the MPRK scheme used is well suited to solve the problem. 
+We can see that the MPRK scheme used is well suited to solve the problem.
 
 
 ## Package versions
diff --git a/docs/src/stratospheric_reaction.md b/docs/src/stratospheric_reaction.md
index a458c9931..4c1d73589 100644
--- a/docs/src/stratospheric_reaction.md
+++ b/docs/src/stratospheric_reaction.md
@@ -54,7 +54,7 @@ In addition, all production and destruction terms not listed have the value zero
 
 ## Solution of the production-destruction system
 
-Now we are ready to define a [`PDSProblem`](@ref) 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/).
+Now we are ready to define a [`PDSProblem`](@ref) and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 
 As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are
 1. an out-of-place implementation with standard (dynamic) matrices and vectors,
@@ -153,7 +153,7 @@ plot(sol_oop,
     widen = true
     )
 ```
-[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) provides the function [`isnonnegative`](@ref) (and also [`isnegative`](@ref)) to check if the solution is actually nonnegative, as expected from an MPRK scheme.
+[PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl) provides the function [`isnonnegative`](@ref) (and also [`isnegative`](@ref)) to check if the solution is actually nonnegative, as expected from an MPRK scheme.
 ```@example stratreac
 isnonnegative(sol_oop)
 ```
@@ -358,7 +358,7 @@ plot(sol_static,
     )
 ```
 
-The above implementation of the stratospheric reaction problem using `StaticArrays` can also be found in the [Example Problems](https://skopecz.github.io/PositiveIntegrators.jl/dev/api_reference/#Example-problems) as [`prob_pds_stratreac`](@ref).
+The above implementation of the stratospheric reaction problem using `StaticArrays` can also be found in the [Example Problems](https://NumericalMathematics.github.io/PositiveIntegrators.jl/dev/api_reference/#Example-problems) as [`prob_pds_stratreac`](@ref).
 
 ### Preservation of linear invariants
 As MPRK schemes do not preserve general linear invariants, especially when applied to non-conservative PDS, we compute and plot the relative errors with respect to both linear invariants to see how well these are preserved.
diff --git a/docs/src/stratospheric_reaction_benchmark.md b/docs/src/stratospheric_reaction_benchmark.md
index 4b746a913..31481430f 100644
--- a/docs/src/stratospheric_reaction_benchmark.md
+++ b/docs/src/stratospheric_reaction_benchmark.md
@@ -1,6 +1,6 @@
 # [Benchmark: Solution of a stratospheric reaction problem](@id benchmark-stratos)
 
-We use the stiff stratospheric reaction problem [`prob_pds_stratreac`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+We use the stiff stratospheric reaction problem [`prob_pds_stratreac`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/NumericalMathematics/PositiveIntegrators.jl).
 
 ```@example stratreac
 using OrdinaryDiffEqFIRK, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqSDIRK
@@ -59,7 +59,7 @@ function stratreac_plot(sols, labels = fill("", length(sols)), sol_ref = nothing
 end
 nothing # hide
 ```
-First, we show approximations of `Rosenbrock23()` using loose tolerances. 
+First, we show approximations of `Rosenbrock23()` using loose tolerances.
 
 ```@example stratreac
 # compute reference solution for plotting
@@ -80,11 +80,11 @@ Although not visible in the plots, the `Rosenbrock23` solution contains negative
 isnonnegative(sol_Ros23)
 ```
 
-Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions. 
+Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions.
 
 ```@example stratreac
 # compute solution with isoutofdomain = isnegative
-sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol, 
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol,
                   isoutofdomain = isnegative); #reject negative solutions
 
 # plot solution
@@ -114,7 +114,7 @@ The remaining poor approximation of the O₂ component could be due to the fact
 
 ## Work-Precision diagrams
 
-In the following we show several work-precision diagrams, which compare different methods with respect to computing times and errors. 
+In the following we show several work-precision diagrams, which compare different methods with respect to computing times and errors.
 First we focus on adaptive methods, afterwards we also show results obtained with fixed time step sizes.
 
 Since the stratospheric reaction problem is stiff, we need to use a suited implicit scheme to compute its reference solution.
@@ -146,14 +146,14 @@ nothing # hide
 
 We also note that MPRK schemes with stricter tolerances, quickly require more than a million time steps, which makes these schemes inefficient in such situations.
 
-First we compare different MPRK schemes. In addition to the default version we also use the schemes with `small_constant = 1e-6`. 
+First we compare different MPRK schemes. In addition to the default version we also use the schemes with `small_constant = 1e-6`.
 
 ```@example stratreac
 # choose methods to compare
 algs = [MPRK22(1.0); MPRK22(1.0, small_constant = 1e-6); SSPMPRK22(0.5, 1.0); SSPMPRK22(0.5, 1.0, small_constant = 1e-6);
         MPRK43I(1.0, 0.5); MPRK43I(1.0, 0.5, small_constant = 1e-6); MPRK43I(0.5, 0.75); MPRK43I(0.5, 0.75, small_constant = 1e-6)
         MPRK43II(0.5); MPRK43II(0.5, small_constant = 1e-6); MPRK43II(2.0 / 3.0); MPRK43II(2.0 / 3.0, small_constant = 1e-6)]
-labels = ["MPRK22(1.0)"; "MPRK22(1.0, sc=1e-6)"; "SSPMPRK22(0.5,1.0)"; "SSPMPRK22(0.5,1.0, sc=1e-6)"; 
+labels = ["MPRK22(1.0)"; "MPRK22(1.0, sc=1e-6)"; "SSPMPRK22(0.5,1.0)"; "SSPMPRK22(0.5,1.0, sc=1e-6)";
           "MPRK43I(1.0,0.5)"; "MPRK43I(1.0,0.5, sc=1e-6)"; "MPRK43I(0.5,0.75)"; "MPRK43I(0.5,0.75, sc=1e-6)"; "MPRK43II(0.5)"; "MPRK43II(0.5, sc=1e-6)"
           "MPRK43II(2.0/3.0)"; "MPRK43II(2.0/3.0, sc=1e-6)"]
 
@@ -161,13 +161,13 @@ labels = ["MPRK22(1.0)"; "MPRK22(1.0, sc=1e-6)"; "SSPMPRK22(0.5,1.0)"; "SSPMPRK2
 wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref; compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,
      color = permutedims([repeat([1],2)..., repeat([2],2)..., repeat([3],4)..., repeat([4],4)...]),
      xlims = (10^-7, 10^0), xticks = 10.0 .^ (-8:1:0),
      ylims = (10^-5, 10^1), yticks = 10.0 .^ (-5:1:1), minorticks = 10)
 ```
 
-We see that using `small_constant = 1e-6` clearly improves the performance of most methods. 
+We see that using `small_constant = 1e-6` clearly improves the performance of most methods.
 For comparisons with other second- and third-order schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) we choose the second-order scheme `MPRK22(1.0, small_constant = 1e-6)` and the third-order scheme `MPRK43I(0.5, 0.75)`.
 To guarantee positive solutions of the [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) methods, we select the solver option `isoutofdomain = isnegative`.
 
@@ -182,11 +182,11 @@ labels2 = ["TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbro
 
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; compute_error)
-work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref; compute_error, 
+work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref; compute_error,
                         isoutofdomain = isnegative)
 
 # plot work-precision diagram
-plot(wp, [labels1; labels2]; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+plot(wp, [labels1; labels2]; title = "Stratospheric reaction benchmark", legend = :bottomleft,
      color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)...]),
      xlims = (10^-8, 10^0), xticks = 10.0 .^ (-8:1:0),
      ylims = (2*10^-4, 5*10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
@@ -203,11 +203,11 @@ labels3 = ["Rodas5P"; "Rodas4P"; "RadauIIA5"]
 
 # compute work-precision data
 wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; compute_error)
-work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref; compute_error, 
+work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref; compute_error,
                         isoutofdomain = isnegative)
 
 # plot work-precision diagram
-plot(wp, [labels1; labels3]; title = "Stratospheric reaction benchmark", legend = :topright,     
+plot(wp, [labels1; labels3]; title = "Stratospheric reaction benchmark", legend = :topright,
      color = permutedims([1, 3, repeat([4], 3)...]),
      xlims = (10^-7, 10^0), xticks = 10.0 .^ (-8:1:0),
      ylims = (2*10^-4, 5*10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
@@ -219,7 +219,7 @@ Again, it can be seen that MPRK methods are only advantageous if low accuracy is
 
 Here we use fixed time step sizes instead of adaptive time stepping.
 We use the functions [`work_precision_fixed`](@ref) and [`work_precision_fixed!`](@ref) to compute the data for the diagrams.
-Please note that these functions set error and computing time to `Inf`, whenever a solution contains negative elements. 
+Please note that these functions set error and computing time to `Inf`, whenever a solution contains negative elements.
 Consequently, such cases are not visible in the work-precision diagrams.
 
 Within the work-precision diagrams we use the following time step sizes.
@@ -258,14 +258,14 @@ Based on the above comparison, we will only consider schemes in which `small_con
 # select schemes
 algs = [MPRK22(1.0); SSPMPRK22(0.5, 1.0); MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0);
         SSPMPRK43()]
-labels = ["MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)"; "MPRK43I(1.0,0.5)"; "MPRK43I(0.5,0.75)";  "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)"; 
+labels = ["MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)"; "MPRK43I(1.0,0.5)"; "MPRK43I(0.5,0.75)";  "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)";
           "SSPMPRK43()"]
 
 # compute work-precision data
 wp = work_precision_fixed(prob, algs, labels, dts, alg_ref; compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,
      color = permutedims([1, 2, repeat([3],2)..., repeat([4],2)..., 5]),
      xlims = (10^-6, 10^2), xticks = 10.0 .^ (-8:1:2),
      ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:1), minorticks = 10)
@@ -277,16 +277,16 @@ For the chosen time step sizes none of the above used standard schemes provides
 
 ```@example stratreac
 # select reference MPRK methods
-algs = [MPRK22(1.0); MPRK43II(0.5); TRBDF2(); Kvaerno3(); KenCarp3(); Rodas3(); ROS2(); ROS3(); Rosenbrock23(); 
+algs = [MPRK22(1.0); MPRK43II(0.5); TRBDF2(); Kvaerno3(); KenCarp3(); Rodas3(); ROS2(); ROS3(); Rosenbrock23();
          Rodas5P(); Rodas4P()]
-labels = ["MPRK22(1.0)"; "MPRK43II(0.5)"; "TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbrock23"; 
+labels = ["MPRK22(1.0)"; "MPRK43II(0.5)"; "TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbrock23";
           "Rodas5P"; "Rodas4P"]
 
 # compute work-precision data
 wp = work_precision_fixed(prob, algs, labels, dts, alg_ref; compute_error)
 
 # plot work-precision diagram
-plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,
      color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 3)...]),
      xlims = (10^-6, 10^1), xticks = 10.0 .^ (-12:2:4),
      ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
diff --git a/test/runtests.jl b/test/runtests.jl
index 64b6c2fa5..8c7e90c8e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -176,7 +176,7 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [
 
 # Function definitions for testset "Linear advection"
 # Functions are defined outside the testset to avoid failing allocation tests,
-# see https://github.com/SKopecz/PositiveIntegrators.jl/pull/89.
+# see https://github.com/NumericalMathematics/PositiveIntegrators.jl/pull/89.
 #
 # in-place syntax for f
 function linear_advection_fd_upwind_f!(du, u, p, t)
@@ -2592,6 +2592,13 @@ end
                 # 2.5 times the mean value. In a loglog plot these
                 # differences won't be significant.
                 @test maximum((v .- m1) ./ m1) < 1.5
+                allowed = 1.5
+                if Sys.isapple()
+                    # On macOS, the time differences are larger
+                    # than on other platforms in CI sometimes.
+                    allowed = 10.0
+                end
+                @test maximum((v .- m1) ./ m1) < allowed
             end
         end
 
@@ -2626,7 +2633,13 @@ end
                     # This test allows computing times that are
                     # 2.5 times the mean value. In a loglog plot these
                     # differences won't be significant.
-                    @test maximum((v .- m1) ./ m1) < 1.5
+                    allowed = 1.5
+                    if Sys.isapple()
+                        # On macOS, the time differences are larger
+                        # than on other platforms in CI sometimes.
+                        allowed = 10.0
+                    end
+                    @test maximum((v .- m1) ./ m1) < allowed
                 end
             end
 
@@ -2657,7 +2670,13 @@ end
                     # This test allows computing times that are
                     # 2.5 times the mean value. In a loglog plot these
                     # differences won't be significant.
-                    @test maximum((v .- m1) ./ m1) < 1.5
+                    allowed = 1.5
+                    if Sys.isapple()
+                        # On macOS, the time differences are larger
+                        # than on other platforms in CI sometimes.
+                        allowed = 10.0
+                    end
+                    @test maximum((v .- m1) ./ m1) < allowed
                 end
             end
         end