From 09133d772cea25a6f95a36919a47e554bc2a8b9d Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Sun, 28 May 2023 08:51:56 -0400
Subject: [PATCH 1/9] overload _alg_autodiff (#238)

---
 Project.toml     | 2 +-
 src/alg_utils.jl | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 46568618a..cd927574e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,7 +38,7 @@ ForwardDiff = "0.10"
 FunctionWrappersWrappers = "0.1.3"
 GaussianDistributions = "0.5"
 Octavian = "0.3.17"
-OrdinaryDiffEq = "6.49.1"
+OrdinaryDiffEq = "6.52"
 PSDMatrices = "0.4.2"
 PrecompileTools = "1"
 RecipesBase = "1"
diff --git a/src/alg_utils.jl b/src/alg_utils.jl
index afea40303..f41aa3ba3 100644
--- a/src/alg_utils.jl
+++ b/src/alg_utils.jl
@@ -3,8 +3,8 @@
 # https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/alg_utils.jl
 ############################################################################################
 
-OrdinaryDiffEq.alg_autodiff(alg::AbstractEK) = true
-OrdinaryDiffEq.alg_autodiff(::EK1{CS,AD}) where {CS,AD} = AD
+OrdinaryDiffEq._alg_autodiff(::AbstractEK) = Val{true}()
+OrdinaryDiffEq._alg_autodiff(::EK1{CS,AD}) where {CS,AD} = Val{AD}()
 OrdinaryDiffEq.alg_difftype(::EK1{CS,AD,DiffType}) where {CS,AD,DiffType} = DiffType
 OrdinaryDiffEq.standardtag(::AbstractEK) = false
 OrdinaryDiffEq.standardtag(::EK1{CS,AD,DiffType,ST}) where {CS,AD,DiffType,ST} = ST

From 7fee244e734a0e383fa797003bf7387ec795af3a Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Sun, 28 May 2023 12:55:59 +0000
Subject: [PATCH 2/9] Set version to 0.11.2

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index cd927574e..52823f72b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ProbNumDiffEq"
 uuid = "bf3e78b0-7d74-48a5-b855-9609533b56a5"
 authors = ["Nathanael Bosch"]
-version = "0.11.1"
+version = "0.11.2"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"

From 1016675d0cb51260d0099fa05c3940726f0d208e Mon Sep 17 00:00:00 2001
From: Nathanael Bosch <nathanael.bosch@uni-tuebingen.de>
Date: Wed, 31 May 2023 13:58:28 +0200
Subject: [PATCH 3/9] Change CI to test LTS+stable+nightly; downstream only on
 stable (#240)

---
 .github/workflows/CI.yml         | 3 ++-
 .github/workflows/Downstream.yml | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 2f5810307..403e0d88f 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -20,7 +20,8 @@ jobs:
       matrix:
         version:
           - '1.6'
-          - '1.8'
+          - '1'
+          - 'nightly'
         os:
           - ubuntu-latest
         arch:
diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml
index 7f9d89155..a4fff7ef4 100644
--- a/.github/workflows/Downstream.yml
+++ b/.github/workflows/Downstream.yml
@@ -15,7 +15,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        julia-version: [1,1.6]
+        julia-version: [1]
         os: [ubuntu-latest]
         package:
           - {user: nathanaelbosch, repo: Fenrir.jl, group: All}

From 646aa33b7298c2184754d976f54a7613fff38283 Mon Sep 17 00:00:00 2001
From: Nathanael Bosch <nathanael.bosch@uni-tuebingen.de>
Date: Thu, 1 Jun 2023 08:49:40 +0200
Subject: [PATCH 4/9] Probabilistic Exponential Integrators (#241)

---
 Project.toml                        |  2 +
 docs/make.jl                        |  1 +
 docs/src/exponential_integrators.md | 95 +++++++++++++++++++++++++++++
 docs/src/index.md                   |  5 +-
 docs/src/solvers.md                 |  8 ++-
 src/ProbNumDiffEq.jl                |  2 +
 src/algorithms.jl                   | 66 ++++++++++++++++++++
 src/caches.jl                       |  5 +-
 src/perform_step.jl                 |  7 +++
 src/priors/common.jl                | 10 +--
 src/priors/ioup.jl                  | 42 +++++++++----
 src/priors/ltisde.jl                | 35 +++++++++++
 test/exponential_integrators.jl     | 41 +++++++++++++
 test/runtests.jl                    |  3 +
 test/stats.jl                       | 33 +++++-----
 15 files changed, 317 insertions(+), 38 deletions(-)
 create mode 100644 docs/src/exponential_integrators.md
 create mode 100644 test/exponential_integrators.jl

diff --git a/Project.toml b/Project.toml
index 52823f72b..663ba3a2c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 GaussianDistributions = "43dcc890-d446-5863-8d1a-14597580bb8d"
@@ -34,6 +35,7 @@ DiffEqBase = "6.122"
 DiffEqDevTools = "2"
 ExponentialUtilities = "1"
 FastBroadcast = "0.2"
+FastGaussQuadrature = "0.5"
 ForwardDiff = "0.10"
 FunctionWrappersWrappers = "0.1.3"
 GaussianDistributions = "0.5"
diff --git a/docs/make.jl b/docs/make.jl
index 60c21fa10..b09afae73 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -11,6 +11,7 @@ makedocs(
             "Getting Started" => "getting_started.md"
             "Second Order ODEs and Energy Preservation" => "dynamical_odes.md"
             "Differential Algebraic Equations" => "dae.md"
+            "Probabilistic Exponential Integrators" => "exponential_integrators.md"
             "Parameter Inference" => "fenrir.md"
         ],
         "Solvers and Options" => [
diff --git a/docs/src/exponential_integrators.md b/docs/src/exponential_integrators.md
new file mode 100644
index 000000000..033e43c04
--- /dev/null
+++ b/docs/src/exponential_integrators.md
@@ -0,0 +1,95 @@
+# Probabilistic Exponential Integrators
+
+Exponential integrators are a class of numerical methods for solving semi-linear ordinary differential equations (ODEs) of the form
+```math
+\begin{aligned}
+\dot{y}(t) &= L y(t) + f(y(t), t), \quad y(0) = y_0,
+\end{aligned}
+```
+where $L$ is a linear operator and $f$ is a nonlinear function.
+In a nutshell, exponential integrators solve the linear part of the ODE exactly, and only approximate the nonlinear part.
+[Probabilistic exponential integrators](https://arxiv.org/abs/2305.14978) are the probabilistic numerics approach to exponential integrators.
+
+## Example
+
+Let's consider a simple semi-linear ODE
+```math
+\begin{aligned}
+\dot{y}(t) &= - y(t) + \sin(y(t)), \quad y(0) = 1.0.
+\end{aligned}
+```
+
+We can solve this ODE reasonably well with the standard [`EK1`](@ref)  and adaptive steps (the default):
+```@example expint
+using ProbNumDiffEq, Plots, LinearAlgebra
+theme(:default; palette=["#4063D8", "#389826", "#9558B2", "#CB3C33"])
+
+f(du, u, p, t) = (@. du = -u + sin(u))
+u0 = [1.0]
+tspan = (0.0, 20.0)
+prob = ODEProblem(f, u0, tspan)
+
+ref = solve(prob, EK1(), abstol=1e-10, reltol=1e-10)
+plot(ref, color=:black, linestyle=:dash, label="Reference")
+```
+
+But for fixed (large) step sizes this ODE is more challenging:
+The explicit [`EK0`](@ref) method oscillates and diverges due to the stiffness of the ODE,
+and the semi-implicit [`EK1`](@ref) method is stable but the solution is not very accurate.
+```@example expint
+STEPSIZE = 4
+DM = FixedDiffusion() # recommended for fixed steps
+
+# we don't smooth the EK0 here to show the oscillations more clearly
+sol0 = solve(prob, EK0(smooth=false, diffusionmodel=DM), adaptive=false, dt=STEPSIZE, dense=false)
+sol1 = solve(prob, EK1(diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
+
+plot(ylims=(0.3, 1.05))
+plot!(ref, color=:black, linestyle=:dash, label="Reference")
+plot!(sol0, denseplot=false, marker=:o, markersize=2, label="EK0", color=1)
+plot!(sol1, denseplot=false, marker=:o, markersize=2, label="EK1", color=2)
+```
+
+_**Probabilistic exponential integrators**_ leverage the semi-linearity of the ODE to compute more accurate solutions for the same fixed step size.
+You can use either the [`ExpEK`](@ref) method and provide the linear part (with the keyword argument `L`),
+or the [`RosenbrockExpEK`](@ref) to automatically linearize along the mean of the numerical solution:
+```@example expint
+sol_exp = solve(prob, ExpEK(L=-1, diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
+sol_ros = solve(prob, RosenbrockExpEK(diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
+
+plot(ylims=(0.3, 1.05))
+plot!(ref, color=:black, linestyle=:dash, label="Reference")
+plot!(sol_exp, denseplot=false, marker=:o, markersize=2, label="ExpEK", color=3)
+plot!(sol_ros, denseplot=false, marker=:o, markersize=2, label="RosenbrockExpEK", color=4)
+```
+
+The solutions are indeed much more accurate than those of the standard `EK1`, for the same fixed step size!
+
+
+## Background: Integrated Ornstein-Uhlenbeck priors
+
+Probabilistic exponential integrators "solve the linear part exactly" by including it into the prior model of the solver.
+Namely, the solver chooses a (q-times) integrated Ornstein-Uhlenbeck prior with rate parameter equal to the linearity.
+The [`ExpEK`](@ref) solver is just a short-hand for an [`EK0`](@ref) with appropriate prior:
+```@repl expint
+ExpEK(order=3, L=-1) == EK0(prior=IOUP(3, -1))
+```
+Similarly, the [`RosenbrockExpEK`](@ref) solver is also just a short-hand:
+```@repl expint
+RosenbrockExpEK(order=3) == EK1(prior=IOUP(3, update_rate_parameter=true))
+```
+
+This means that you can also construct other probabilistic exponential integrators by hand!
+In this example the `EK1` with `IOUP` prior with rate parameter `-1` performs extremely well:
+```@example expint
+sol_expek1 = solve(prob, EK1(prior=IOUP(3, -1), diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
+
+plot(ylims=(0.3, 1.05))
+plot!(ref, color=:black, linestyle=:dash, label="Reference")
+plot!(sol_expek1, denseplot=false, marker=:o, markersize=2, label="EK1 + IOUP")
+```
+
+
+### References
+
+[1] N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023) ([link](https://arxiv.org/abs/2305.14978))
diff --git a/docs/src/index.md b/docs/src/index.md
index 1012d4a52..c551b277a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -45,8 +45,9 @@ For a quick introduction check out the "[Solving ODEs with Probabilistic Numeric
 
 ## [References](@id references)
 
+  - N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023)
   - N. Bosch, F. Tronarp, P. Hennig: **Pick-and-Mix Information Operators for Probabilistic ODE Solvers** (2022)
-  - N. Krämer, N. Bosch, J. Schmidt, P. Hennig: **Probabilistic ODE Solutions in Millions of Dimensions** (2021)
+  - N. Krämer, N. Bosch, J. Schmidt, P. Hennig: **Probabilistic ODE Solutions in Millions of Dimensions** (2022)
   - N. Bosch, P. Hennig, F. Tronarp: **Calibrated Adaptive Probabilistic ODE Solvers** (2021)
   - F. Tronarp, S. Särkkä, and P. Hennig: **Bayesian ODE Solvers: The Maximum A Posteriori Estimate** (2021)
   - N. Krämer, P. Hennig: **Stable Implementation of Probabilistic ODE Solvers** (2020)
@@ -54,4 +55,4 @@ For a quick introduction check out the "[Solving ODEs with Probabilistic Numeric
   - F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig: **Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective** (2019)
   - M. Schober, S. Särkkä, and P. Hennig: **A Probabilistic Model for the Numerical Solution of Initial Value Problems** (2018)
 
-A much more detailed list of references, not only on ODE filters but on probabilistic numerics in general, can be found on the [probabilistic-numerics.org homepage](https://www.probabilistic-numerics.org/research/general/).
+More references on ODE filters and on probabilistic numerics in general can be found on the [probabilistic-numerics.org homepage](https://www.probabilistic-numerics.org/research/general/).
diff --git a/docs/src/solvers.md b/docs/src/solvers.md
index 5620a593b..62eb5844c 100644
--- a/docs/src/solvers.md
+++ b/docs/src/solvers.md
@@ -7,14 +7,18 @@ All solvers are compatible with DAEs in mass-matrix ODE form.
 They also specialize on second-order ODEs: If the problem is of type [`SecondOrderODEProblem`](https://docs.sciml.ai/DiffEqDocs/stable/types/dynamical_types/#SciMLBase.SecondOrderODEProblem), it solves the second-order problem directly; this is more efficient than solving the transformed first-order problem and provides more meaningful posteriors
 [[1]](@ref solversrefs).
 
-
 ## API
-
 ```@docs
 EK1
 EK0
 ```
 
+### Probabilistic Exponential Integrators
+```@docs
+ExpEK
+RosenbrockExpEK
+```
+
 ## [References](@id solversrefs)
 
 [1] N. Bosch, F. Tronarp, P. Hennig: **Pick-and-Mix Information Operators for Probabilistic ODE Solvers** (2022) ([link](https://proceedings.mlr.press/v151/bosch22a.html))
diff --git a/src/ProbNumDiffEq.jl b/src/ProbNumDiffEq.jl
index dda8a02fe..109ab007f 100644
--- a/src/ProbNumDiffEq.jl
+++ b/src/ProbNumDiffEq.jl
@@ -25,6 +25,7 @@ using RecursiveArrayTools
 using ForwardDiff
 using ExponentialUtilities
 using Octavian
+using FastGaussQuadrature
 
 @reexport using GaussianDistributions
 using GaussianDistributions: logpdf
@@ -57,6 +58,7 @@ export TaylorModeInit, ClassicSolverInit
 
 include("algorithms.jl")
 export EK0, EK1
+export ExpEK, RosenbrockExpEK
 
 include("alg_utils.jl")
 include("caches.jl")
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 79341c0d3..eb53e6251 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -111,6 +111,72 @@ EK1(;
         initialization,
     )
 
+"""
+    ExpEK(; L, order=3, kwargs...)
+
+Probabilistic exponential integrator
+
+Probabilistic exponential integrators are a class of integrators for semi-linear stiff ODEs
+that provide improved stability by essentially solving the linear part of the ODE exactly.
+In probabilistic numerics, this amounts to including the linear part into the prior model
+of the solver.
+
+`ExpEK` is therefore just a short-hand for [`EK0`](@ref) with [`IOUP`](@ref) prior:
+```julia
+ExpEK(; order=3, L, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)
+```
+
+See also [`RosenbrockExpEK`](@ref), [`EK0`](@ref), [`EK1`](@ref).
+
+# Arguments
+See [`EK0`](@ref) for available keyword arguments.
+
+# Examples
+```julia-repl
+julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
+julia> solve(prob, ExpEK(L=-1))
+```
+
+
+# Reference
+[1] N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023) ([link](https://arxiv.org/abs/2305.14978))
+"""
+ExpEK(; L, order=3, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)
+
+"""
+    RosenbrockExpEK(; order=3, kwargs...)
+
+**Probabilistic Rosenbrock-type exponential integrator**
+
+A probabilistic exponential integrator similar to [`ExpEK`](@ref), but with automatic
+linearization along the mean numerical solution. This brings the advantage that the
+linearity does not need to be specified manually, and the more accurate local linearization
+can sometimes also improve stability; but since the "prior" is adjusted at each step the
+probabilistic interpretation becomes more complicated.
+
+`RosenbrockExpEK` is just a short-hand for [`EK1`](@ref) with appropriete [`IOUP`](@ref)
+prior:
+```julia
+RosenbrockExpEK(; order=3, kwargs...) = EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)
+```
+
+See also [`ExpEK`](@ref), [`EK0`](@ref), [`EK1`](@ref).
+
+# Arguments
+See [`EK1`](@ref) for available keyword arguments.
+
+# Examples
+```julia-repl
+julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
+julia> solve(prob, RosenbrockExpEK())
+```
+
+# Reference
+[1] N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023) ([link](https://arxiv.org/abs/2305.14978))
+"""
+RosenbrockExpEK(; order=3, kwargs...) =
+    EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)
+
 function DiffEqBase.remake(thing::EK1{CS,AD,DT,ST,CJ}; kwargs...) where {CS,AD,DT,ST,CJ}
     T = SciMLBase.remaker_of(thing)
     T(;
diff --git a/src/caches.jl b/src/caches.jl
index 498b23f16..9f0143693 100644
--- a/src/caches.jl
+++ b/src/caches.jl
@@ -105,8 +105,11 @@ function OrdinaryDiffEq.alg_cache(
     # Prior dynamics
     prior = if alg.prior isa IWP
         IWP{uElType}(d, alg.prior.num_derivatives)
+    elseif alg.prior isa IOUP && ismissing(alg.prior.rate_parameter)
+        r = zeros(uElType, d, d)
+        IOUP{uElType}(d, q, r, alg.prior.update_rate_parameter)
     elseif alg.prior isa IOUP
-        IOUP{uElType}(d, q, alg.prior.rate_parameter)
+        IOUP{uElType}(d, q, alg.prior.rate_parameter, alg.prior.update_rate_parameter)
     elseif alg.prior isa Matern
         Matern{uElType}(d, q, alg.prior.lengthscale)
     else
diff --git a/src/perform_step.jl b/src/perform_step.jl
index cb695c956..2e56af026 100644
--- a/src/perform_step.jl
+++ b/src/perform_step.jl
@@ -24,6 +24,8 @@ function make_new_transitions(integ, cache, repeat_step)::Bool
     # Similar to OrdinaryDiffEq.do_newJ
     if integ.iter <= 1
         return true
+    elseif cache.prior isa IOUP && cache.prior.update_rate_parameter
+        return true
     elseif repeat_step
         return false
     elseif integ.dt == cache.dt_last
@@ -60,6 +62,11 @@ function OrdinaryDiffEq.perform_step!(integ, cache::EKCache, repeat_step=false)
     tnew = t + dt
 
     if make_new_transitions(integ, cache, repeat_step)
+        # Rosenbrock-style update of the IOUP rate parameter
+        if cache.prior isa IOUP && cache.prior.update_rate_parameter
+            OrdinaryDiffEq.calc_J!(cache.prior.rate_parameter, integ, cache, false)
+        end
+
         make_transition_matrices!(cache, cache.prior, dt)
     end
 
diff --git a/src/priors/common.jl b/src/priors/common.jl
index 6d91fa01a..adb142e54 100644
--- a/src/priors/common.jl
+++ b/src/priors/common.jl
@@ -32,10 +32,12 @@ See also: [`make_transition_matrices`](@ref).
 [1] N. Krämer, P. Hennig: **Stable Implementation of Probabilistic ODE Solvers** (2020)
 """
 function initialize_transition_matrices(p::AbstractODEFilterPrior{T}, dt) where {T}
-    Ah, Qh = discretize(p, dt)
+    d, q = p.wiener_process_dimension, p.num_derivatives
+    D = d * (q + 1)
+    Ah, Qh = zeros(T, D, D), PSDMatrix(zeros(T, D, D))
     P, PI = initialize_preconditioner(p, dt)
-    A = P * Ah * PI
-    Q = X_A_Xt(Qh, P)
+    A = copy(Ah)
+    Q = copy(Qh)
     return A, Q, Ah, Qh, P, PI
 end
 
@@ -71,7 +73,7 @@ function make_transition_matrices!(cache, prior::AbstractODEFilterPrior, dt)
     _Ah, _Qh = discretize(cache.prior, dt)
     copy!(Ah, _Ah)
     copy!(Qh, _Qh)
-    A .= P.diag .* Ah .* PI.diag'
+    @.. A = P.diag * Ah * PI.diag'
     fast_X_A_Xt!(Q, Qh, P)
 end
 
diff --git a/src/priors/ioup.jl b/src/priors/ioup.jl
index 4a4628dbc..0840138d7 100644
--- a/src/priors/ioup.jl
+++ b/src/priors/ioup.jl
@@ -28,16 +28,37 @@ struct IOUP{elType,dimType,R} <: AbstractODEFilterPrior{elType}
     wiener_process_dimension::dimType
     num_derivatives::Int
     rate_parameter::R
+    update_rate_parameter::Bool
 end
-IOUP(num_derivatives, rate_parameter) =
-    IOUP(missing, num_derivatives, rate_parameter)
-IOUP(wiener_process_dimension, num_derivatives, rate_parameter) =
-    IOUP{typeof(1.0)}(wiener_process_dimension, num_derivatives, rate_parameter)
-IOUP{T}(wiener_process_dimension, num_derivatives, rate_parameter) where {T} =
+IOUP(num_derivatives; update_rate_parameter) = begin
+    @assert update_rate_parameter
+    IOUP(missing, num_derivatives, missing, update_rate_parameter)
+end
+IOUP(num_derivatives, rate_parameter; update_rate_parameter=false) =
+    IOUP(missing, num_derivatives, rate_parameter, update_rate_parameter)
+IOUP(
+    wiener_process_dimension,
+    num_derivatives,
+    rate_parameter,
+    update_rate_parameter=false,
+) =
+    IOUP{typeof(1.0)}(
+        wiener_process_dimension,
+        num_derivatives,
+        rate_parameter,
+        update_rate_parameter,
+    )
+IOUP{T}(
+    wiener_process_dimension,
+    num_derivatives,
+    rate_parameter,
+    update_rate_parameter=false,
+) where {T} =
     IOUP{T,typeof(wiener_process_dimension),typeof(rate_parameter)}(
         wiener_process_dimension,
         num_derivatives,
         rate_parameter,
+        update_rate_parameter,
     )
 
 function to_1d_sde(p::IOUP)
@@ -93,21 +114,16 @@ end
 function discretize(p::IOUP, dt::Real)
     r = p.rate_parameter
     A, Q = if p.rate_parameter isa Number
-        A_breve, Q_breve = discretize(to_1d_sde(p), dt)
-        d = p.wiener_process_dimension
-        # QR_breve = cholesky!(Symmetric(Q_breve)).L'
-        E = eigen(Symmetric(Q_breve))
-        QR_breve = Diagonal(sqrt.(max.(E.values, 0))) * E.vectors'
+        A_breve, QR_breve = discretize_sqrt(to_1d_sde(p), dt)
 
+        d = p.wiener_process_dimension
         A = kron(I(d), A_breve)
         QR = kron(I(d), QR_breve)
         Q = PSDMatrix(QR)
         A, Q
     else
         @assert r isa AbstractVector || r isa AbstractMatrix
-        A, Q = discretize(to_sde(p), dt)
-        E = eigen(Symmetric(Q))
-        QR = Diagonal(sqrt.(max.(E.values, 0))) * E.vectors'
+        A, QR = discretize_sqrt(to_sde(p), dt)
         Q = PSDMatrix(QR)
         A, Q
     end
diff --git a/src/priors/ltisde.jl b/src/priors/ltisde.jl
index 4287bc000..30ba4efda 100644
--- a/src/priors/ltisde.jl
+++ b/src/priors/ltisde.jl
@@ -36,3 +36,38 @@ function matrix_fraction_decomposition(
     Q = Mexp[1:d, d+1:end] * A'
     return A, Q
 end
+
+function discretize_sqrt(sde::LTISDE, dt::Real)
+    F, L = drift(sde), dispersion(sde)
+
+    D = size(F, 1)
+    d = size(L, 2)
+    N = Int(D / d)
+    R = similar(F, N * d, D)
+    method = ExpMethodHigham2005()
+    expcache = ExponentialUtilities.alloc_mem(F, method)
+
+    Ah = exponential!(dt * F, method, expcache)
+
+    chol_integrand(τ) = begin
+        E = exponential!((dt - τ) * F', method, expcache)
+        L'E
+    end
+    nodes, weights = gausslegendre(N)
+    b, a = dt, 0
+    @. nodes = (b - a) / 2 * nodes + (a + b) / 2
+    @. weights = (b - a) / 2 * weights
+    @simd ivdep for i in 1:N
+        R[(i-1)*d+1:i*d, 1:D] .= sqrt(weights[i]) .* chol_integrand(nodes[i])
+    end
+
+    M = R'R |> Symmetric
+    chol = cholesky!(M, check=false)
+    Qh_R = if issuccess(chol)
+        chol.U |> Matrix
+    else
+        qr!(R).R |> Matrix
+    end
+
+    return Ah, Qh_R
+end
diff --git a/test/exponential_integrators.jl b/test/exponential_integrators.jl
new file mode 100644
index 000000000..0c7dcfae6
--- /dev/null
+++ b/test/exponential_integrators.jl
@@ -0,0 +1,41 @@
+using ProbNumDiffEq
+using LinearAlgebra
+using Random
+using Test
+
+@testset "Linear ODE" begin
+    f!(du, u, p, t) = @. du = p * u
+    u0 = [1.0]
+    tspan = (0.0, 10.0)
+    p = -1
+    prob = ODEProblem(f!, u0, tspan, p)
+
+    u(t) = u0 * exp(p * t)
+    uend = u(tspan[2])
+
+    sol0 = solve(prob, EK0(order=3))
+    sol1 = solve(prob, EK1(order=3))
+    solexp = solve(prob, ExpEK(L=p, order=3))
+    solros = solve(prob, RosenbrockExpEK(order=3))
+
+    err0 = norm(uend - sol0[end])
+    @test err0 < 1e-7
+    err1 = norm(uend - sol1[end])
+    @test err1 < 1e-9
+    errexp = norm(uend - solexp[end])
+    @test errexp < 1e-10
+    errros = norm(uend - solros[end])
+    @test errros < 1e-13
+
+    @test errros < errexp < err1 < err0
+
+    @test length(solexp) < length(sol0)
+    @test solexp.stats.nf < sol0.stats.nf
+    @test length(solexp) < length(sol1)
+    @test solexp.stats.nf < sol1.stats.nf
+
+    @test length(solros) < length(sol0)
+    @test solros.stats.nf < sol0.stats.nf
+    @test length(solros) < length(sol1)
+    @test solros.stats.nf < sol1.stats.nf
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index a64bf9368..69b94d53c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -58,6 +58,9 @@ const GROUP = get(ENV, "GROUP", "All")
             @timedtestset "IOUP" begin
                 include("ioup.jl")
             end
+            @timedtestset "Exponential Integrators" begin
+                include("exponential_integrators.jl")
+            end
         end
     end
 
diff --git a/test/stats.jl b/test/stats.jl
index 8f206ec60..e9ec43800 100644
--- a/test/stats.jl
+++ b/test/stats.jl
@@ -4,13 +4,15 @@ using LinearAlgebra
 
 const q = 3
 
-@testset "stats.nf testing $alg" for init in (TaylorModeInit(), ClassicSolverInit()),
-    alg in (
-        EK0(order=q, smooth=false, initialization=init),
-        EK1(order=q, smooth=false, initialization=init),
-        EK1(order=q, smooth=false, initialization=init, autodiff=false),
-    )
-
+@testset "stats.nf testing $alg" for alg in (
+    EK0(prior=IWP(q), smooth=false),
+    EK0(prior=IWP(q), smooth=false, initialization=ClassicSolverInit()),
+    EK1(prior=IWP(q), smooth=false),
+    EK1(prior=IWP(q), smooth=false, autodiff=false),
+    EK1(prior=IOUP(q, -1), smooth=false),
+    EK1(prior=IOUP(q, update_rate_parameter=true), smooth=false),
+    EK1(prior=IOUP(q, update_rate_parameter=true), smooth=true),
+)
     f_counter = [0]
     function f(du, u, p, t; f_counter=f_counter)
         f_counter .+= 1
@@ -21,18 +23,17 @@ const q = 3
     p = [-0.1]
     tspan = (0.0, 1.0)
     prob = ODEProblem(f, u0, tspan, p)
-    sol = solve(prob, alg, save_everystep=false, dense=false)
+    sol = solve(prob, alg, save_everystep=alg.smooth, dense=alg.smooth)
     @test sol.stats.nf == f_counter[1]
 end
 
-@testset "SecondOrderODEProblem stats.nf testing $alg" for init in (TaylorModeInit(),),
+@testset "SecondOrderODEProblem stats.nf testing $alg" for alg in (
+    EK0(prior=IWP(q), smooth=false),
+    # EK0(prior=IWP(q), smooth=false, initialization=ClassicSolverInit()),
     # ClassicSolverInit does not work for second order ODEs right now
-    alg in (
-        EK0(order=q, smooth=false, initialization=init),
-        EK1(order=q, smooth=false, initialization=init),
-        EK1(order=q, smooth=false, initialization=init, autodiff=false),
-    )
-
+    EK1(prior=IWP(q), smooth=false),
+    EK1(prior=IWP(q), smooth=false, autodiff=false),
+)
     f_counter = [0]
     function f(ddu, du, u, p, t; f_counter=f_counter)
         f_counter .+= 1
@@ -44,6 +45,6 @@ end
     p = [-1.0]
     tspan = (0.0, 1.0)
     prob = SecondOrderODEProblem(f, du0, u0, tspan, p)
-    sol = solve(prob, alg, save_everystep=false, dense=false)
+    sol = solve(prob, alg, save_everystep=alg.smooth, dense=alg.smooth)
     @test sol.stats.nf == f_counter[1]
 end

From cc3dfd26d6f78c9aacbed01e3e6981cd36085171 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Thu, 1 Jun 2023 06:54:53 +0000
Subject: [PATCH 5/9] Set version to 0.12.0

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 663ba3a2c..7cb2757e4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ProbNumDiffEq"
 uuid = "bf3e78b0-7d74-48a5-b855-9609533b56a5"
 authors = ["Nathanael Bosch"]
-version = "0.11.2"
+version = "0.12.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"

From d2a9bf9a9f08a982efb0be9f0f19958a3d9ab45d Mon Sep 17 00:00:00 2001
From: Nathanael Bosch <nathanael.bosch@uni-tuebingen.de>
Date: Fri, 2 Jun 2023 13:08:14 +0200
Subject: [PATCH 6/9] Documentation updates (#242)

---
 docs/make.jl                                  | 14 ++--
 docs/src/index.md                             | 81 ++++++++++++-------
 docs/src/initialization.md                    |  4 +-
 docs/src/{ => tutorials}/dae.md               |  0
 docs/src/{ => tutorials}/dynamical_odes.md    |  0
 .../exponential_integrators.md                |  0
 docs/src/{ => tutorials}/fenrir.md            |  0
 docs/src/{ => tutorials}/getting_started.md   |  0
 8 files changed, 60 insertions(+), 39 deletions(-)
 rename docs/src/{ => tutorials}/dae.md (100%)
 rename docs/src/{ => tutorials}/dynamical_odes.md (100%)
 rename docs/src/{ => tutorials}/exponential_integrators.md (100%)
 rename docs/src/{ => tutorials}/fenrir.md (100%)
 rename docs/src/{ => tutorials}/getting_started.md (100%)

diff --git a/docs/make.jl b/docs/make.jl
index b09afae73..a2c0680e8 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -8,17 +8,17 @@ makedocs(
     pages=[
         "Home" => "index.md",
         "Tutorials" => [
-            "Getting Started" => "getting_started.md"
-            "Second Order ODEs and Energy Preservation" => "dynamical_odes.md"
-            "Differential Algebraic Equations" => "dae.md"
-            "Probabilistic Exponential Integrators" => "exponential_integrators.md"
-            "Parameter Inference" => "fenrir.md"
+            "Getting Started" => "tutorials/getting_started.md"
+            "Second Order ODEs and Energy Preservation" => "tutorials/dynamical_odes.md"
+            "Differential Algebraic Equations" => "tutorials/dae.md"
+            "Probabilistic Exponential Integrators" => "tutorials/exponential_integrators.md"
+            "Parameter Inference" => "tutorials/fenrir.md"
         ],
         "Solvers and Options" => [
             "solvers.md",
-            "diffusions.md",
-            "initialization.md",
             "priors.md",
+            "initialization.md",
+            "diffusions.md",
         ],
         "Benchmarks" => [
             "Multi-Language Wrapper Benchmark" => "benchmarks/multi-language-wrappers.md",
diff --git a/docs/src/index.md b/docs/src/index.md
index c551b277a..392a7762f 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -8,20 +8,13 @@ The implemented _ODE filters_ solve differential equations via Bayesian filterin
 
 For a short intro video, check out our [poster presentation at JuliaCon2021](https://www.youtube.com/watch?v=EMFl6ytP3iQ).
 
----
-
-__For more probabilistic numerics check out the [ProbNum](https://probnum.readthedocs.io/en/latest/) Python package.__
-It implements probabilistic ODE solvers, but also probabilistic linear solvers, Bayesian quadrature, and many filtering and smoothing implementations.
-
----
-
 ## Installation
 
 Run Julia, enter `]` to bring up Julia's package manager, and add the ProbNumDiffEq.jl package:
 
 ```
 julia> ]
-(v1.8) pkg> add ProbNumDiffEq
+(v1.9) pkg> add ProbNumDiffEq
 ```
 
 ## Getting Started
@@ -30,29 +23,57 @@ For a quick introduction check out the "[Solving ODEs with Probabilistic Numeric
 
 ## Features
 
-  - Two extended Kalman filtering-based probabilistic solvers: the explicit [`EK0`](@ref) and semi-implicit [`EK1`](@ref).
-  - Adaptive step-size selection (by default with PI control)
-  - On-line uncertainty calibration, for multiple different measurement models
-  - Dense output
-  - Sampling from the solution
-  - Callback support
-  - Convenient plotting through a Plots.jl recipe
-  - Automatic differentiation via ForwardDiff.jl
-  - Supports arbitrary precision numbers via BigFloats.jl
-  - Specialized solvers for second-order ODEs (demo will be added)
-  - Compatible with DAEs in mass-matrix ODE form (demo will be added)
+- Two extended Kalman filtering-based probabilistic solvers: the explicit [`EK0`](@ref) and semi-implicit [`EK1`](@ref).
+- Adaptive step-size selection with PI control;
+  fully compatible with [DifferentialEquations.jl's timestepping options](https://docs.sciml.ai/DiffEqDocs/stable/extras/timestepping/)
+- Online uncertainty calibration for multiple different diffusion models (see "[Diffusion models and calibration](@ref)")
+- [Dense output](@ref)
+- Sampling from the solution
+- Callback support
+- Convenient plotting through a [Plots.jl](https://docs.juliaplots.org/latest/) recipe
+- Automatic differentiation via [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
+- Arbitrary precision via Julia's built-in [arbitrary precision arithmetic](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic)
+- Specialized solvers for second-order ODEs (see [Second Order ODEs and Energy Preservation](@ref))
+- Compatible with DAEs in mass-matrix ODE form (see [Solving DAEs with Probabilistic Numerics](@ref))
 
 
 ## [References](@id references)
 
-  - N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023)
-  - N. Bosch, F. Tronarp, P. Hennig: **Pick-and-Mix Information Operators for Probabilistic ODE Solvers** (2022)
-  - N. Krämer, N. Bosch, J. Schmidt, P. Hennig: **Probabilistic ODE Solutions in Millions of Dimensions** (2022)
-  - N. Bosch, P. Hennig, F. Tronarp: **Calibrated Adaptive Probabilistic ODE Solvers** (2021)
-  - F. Tronarp, S. Särkkä, and P. Hennig: **Bayesian ODE Solvers: The Maximum A Posteriori Estimate** (2021)
-  - N. Krämer, P. Hennig: **Stable Implementation of Probabilistic ODE Solvers** (2020)
-  - H. Kersting, T. J. Sullivan, and P. Hennig: **Convergence Rates of Gaussian Ode Filters** (2020)
-  - F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig: **Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective** (2019)
-  - M. Schober, S. Särkkä, and P. Hennig: **A Probabilistic Model for the Numerical Solution of Initial Value Problems** (2018)
-
-More references on ODE filters and on probabilistic numerics in general can be found on the [probabilistic-numerics.org homepage](https://www.probabilistic-numerics.org/research/general/).
+The main references _for this package_ include:
+
+- N. Bosch, P. Hennig, F. Tronarp:
+  [**Probabilistic Exponential Integrators**](https://arxiv.org/abs/2305.14978)
+  (2023)
+- N. Bosch, F. Tronarp, P. Hennig:
+  [**Pick-and-Mix Information Operators for Probabilistic ODE Solvers**](https://proceedings.mlr.press/v151/bosch22a.html)
+  (2022)
+- N. Krämer, N. Bosch, J. Schmidt, P. Hennig:
+  [**Probabilistic ODE Solutions in Millions of Dimensions**](https://proceedings.mlr.press/v162/kramer22b.html)
+  (2022)
+- N. Bosch, P. Hennig, F. Tronarp:
+  [**Calibrated Adaptive Probabilistic ODE Solvers**](http://proceedings.mlr.press/v130/bosch21a.html)
+  (2021)
+- F. Tronarp, S. Särkkä, and P. Hennig:
+  [**Bayesian ODE Solvers: The Maximum A Posteriori Estimate**](https://link.springer.com/article/10.1007/s11222-021-09993-7)
+  (2021)
+- N. Krämer, P. Hennig:
+  [**Stable Implementation of Probabilistic ODE Solvers**](https://arxiv.org/abs/2012.10106v1)
+  (2020)
+- H. Kersting, T. J. Sullivan, and P. Hennig:
+  [**Convergence Rates of Gaussian Ode Filters**](https://link.springer.com/article/10.1007/s11222-020-09972-4)
+  (2020)
+- F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig:
+  [**Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective**](https://link.springer.com/article/10.1007/s11222-019-09900-1)
+  (2019)
+- M. Schober, S. Särkkä, and P. Hennig:
+  [**A Probabilistic Model for the Numerical Solution of Initial Value Problems**](https://link.springer.com/article/10.1007/s11222-017-9798-7)
+  (2018)
+
+More references on ODE filters and on probabilistic numerics in general can be found on [probabilistic-numerics.org ](https://www.probabilistic-numerics.org/research/general/).
+
+
+## Related packages
+
+- [probdiffeq](https://pnkraemer.github.io/probdiffeq/): Fast and feature-rich filtering-based probabilistic ODE solvers in JAX.
+- [ProbNum](https://probnum.readthedocs.io/en/latest/): Probabilistic numerics in Python. It has not only probabilistic ODE solvers, but also probabilistic linear solvers, Bayesian quadrature, and many filtering and smoothing implementations.
+- [Fenrir.jl](https://github.com/nathanaelbosch/Fenrir.jl): Parameter-inference in ODEs with probabilistic ODE solvers. This package builds on ProbNumDiffEq.jl to provide a negative marginal log-likelihood function, which can then be used with an optimizer or with MCMC for parameter inference.
diff --git a/docs/src/initialization.md b/docs/src/initialization.md
index 0e2741b7e..3ed5f2bc7 100644
--- a/docs/src/initialization.md
+++ b/docs/src/initialization.md
@@ -48,8 +48,8 @@ Y^{(0)}(0) &= y_0, \\
 Y^{(1)}(0) &= f(y_0, 0).
 \end{aligned}
 ```
-And it turns out that we can also compute higher-order derivatives of ``y`` with the chain rule,
-which we can use to initialize ``Y^{(i)}(0)``.
+It turns out that we can also compute higher-order derivatives of ``y`` with the chain rule,
+and then use these to better initialize ``Y^{(i)}(0)``.
 This, done efficiently with Taylor-mode autodiff by using
 [TaylorIntegration.jl](https://perezhz.github.io/TaylorIntegration.jl/latest/),
 is what [`TaylorModeInit`](@ref) does.
diff --git a/docs/src/dae.md b/docs/src/tutorials/dae.md
similarity index 100%
rename from docs/src/dae.md
rename to docs/src/tutorials/dae.md
diff --git a/docs/src/dynamical_odes.md b/docs/src/tutorials/dynamical_odes.md
similarity index 100%
rename from docs/src/dynamical_odes.md
rename to docs/src/tutorials/dynamical_odes.md
diff --git a/docs/src/exponential_integrators.md b/docs/src/tutorials/exponential_integrators.md
similarity index 100%
rename from docs/src/exponential_integrators.md
rename to docs/src/tutorials/exponential_integrators.md
diff --git a/docs/src/fenrir.md b/docs/src/tutorials/fenrir.md
similarity index 100%
rename from docs/src/fenrir.md
rename to docs/src/tutorials/fenrir.md
diff --git a/docs/src/getting_started.md b/docs/src/tutorials/getting_started.md
similarity index 100%
rename from docs/src/getting_started.md
rename to docs/src/tutorials/getting_started.md

From 9896a407834c7134140f107f8bfbbb827b318cbe Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Tue, 4 Jul 2023 14:17:57 +0200
Subject: [PATCH 7/9] CompatHelper: bump compat for TaylorIntegration to 0.14,
 (keep existing compat) (#244)

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 7cb2757e4..92e2366ba 100644
--- a/Project.toml
+++ b/Project.toml
@@ -51,7 +51,7 @@ SimpleUnPack = "1"
 SpecialMatrices = "3"
 StaticArrayInterface = "1.3"
 StructArrays = "0.4, 0.5, 0.6"
-TaylorIntegration = "0.8, 0.9, 0.10, 0.11"
+TaylorIntegration = "0.8, 0.9, 0.10, 0.11, 0.14"
 TaylorSeries = "0.10, 0.11, 0.12, 0.13, 0.14, 0.15"
 ToeplitzMatrices = "0.7, 0.8"
 julia = "1.6"

From 221cd1f071ee46fccebc616d1b8dda70feb942a4 Mon Sep 17 00:00:00 2001
From: Nathanael Bosch <nathanael.bosch@uni-tuebingen.de>
Date: Tue, 1 Aug 2023 09:23:58 +0000
Subject: [PATCH 8/9] Use DocumenterCitations.jl to manage references (#245)

---
 docs/Project.toml                             |   2 +
 docs/make.jl                                  |  29 ++-
 docs/src/assets/citations.css                 |  17 ++
 docs/src/diffusions.md                        |   9 +-
 docs/src/index.md                             |  35 ----
 docs/src/initialization.md                    |  10 +-
 docs/src/references.md                        |   5 +
 docs/src/refs.bib                             | 172 ++++++++++++++++++
 docs/src/solvers.md                           |   8 +-
 docs/src/tutorials/dae.md                     |   7 +-
 docs/src/tutorials/dynamical_odes.md          |  14 +-
 docs/src/tutorials/exponential_integrators.md |   9 +-
 docs/src/tutorials/fenrir.md                  |  11 +-
 expdev.jl                                     |  24 +++
 src/algorithms.jl                             |   4 +-
 src/diffusions.jl                             |  10 +-
 src/initialization/common.jl                  |   6 +-
 17 files changed, 312 insertions(+), 60 deletions(-)
 create mode 100644 docs/src/assets/citations.css
 create mode 100644 docs/src/references.md
 create mode 100644 docs/src/refs.bib
 create mode 100644 expdev.jl

diff --git a/docs/Project.toml b/docs/Project.toml
index f1aaf6d6e..1eeef1036 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,7 +1,9 @@
 [deps]
+Bibliography = "f1be7e48-bf82-45af-a471-ae754a193061"
 DiffEqUncertainty = "ef61062a-5684-51dc-bb67-a0fcdec5c97d"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Fenrir = "e9b4b195-f5cd-427c-8076-5358c553c37f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
diff --git a/docs/make.jl b/docs/make.jl
index a2c0680e8..3bf7371ee 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,9 +1,35 @@
 using Documenter
 using ProbNumDiffEq
 
+using DocumenterCitations, Bibliography
+
+# DocumenterCitations.bib_sorting(::Val{:numeric}) = :nyt
+# function DocumenterCitations.format_bibliography_label(
+#     ::Val{:numeric},
+#     entry,
+#     citations::DocumenterCitations.OrderedDict{String,Int64},
+# )
+#     key = entry.id
+#     sorted_bibtex_keys = citations |> keys |> collect |> sort
+#     i = findfirst(x -> x == key, sorted_bibtex_keys)
+#     @info "format_bibliography_label" entry.id citations sorted_bibtex_keys i entry.date
+#     return "[$i]"
+# end
+# Bibliography.sorting_rules[:nyt] = [:authors; :date; :title]
+
+bib = CitationBibliography(
+    joinpath(@__DIR__, "src", "refs.bib"),
+    style=:numeric,
+    # style=:authoryear,
+)
+sort_bibliography!(bib.entries, :nyt)  # name-year-title
+
 makedocs(
+    bib,
     sitename="ProbNumDiffEq.jl",
-    format=Documenter.HTML(),
+    format=Documenter.HTML(
+        assets=String["assets/citations.css"],
+    ),
     modules=[ProbNumDiffEq],
     pages=[
         "Home" => "index.md",
@@ -30,6 +56,7 @@ makedocs(
             "Filtering and Smoothing" => "filtering.md"
             "Implementation via OrdinaryDiffEq.jl" => "implementation.md"
         ],
+        "References" => "references.md",
     ],
 )
 
diff --git a/docs/src/assets/citations.css b/docs/src/assets/citations.css
new file mode 100644
index 000000000..40a19a310
--- /dev/null
+++ b/docs/src/assets/citations.css
@@ -0,0 +1,17 @@
+.citation dl {
+    display: grid;
+    grid-template-columns: max-content auto; }
+.citation dt {
+    grid-column-start: 1; }
+.citation dd {
+    grid-column-start: 2;
+    margin-bottom: 0.75em; }
+.citation ul {
+    padding: 0 0 2.25em 0;
+    margin: 0;
+    list-style: none;}
+.citation ul li {
+    text-indent: -2.25em;
+    margin: 0.33em 0.5em 0.5em 2.25em;}
+.citation ol li {
+    padding-left:0.75em;}
diff --git a/docs/src/diffusions.md b/docs/src/diffusions.md
index 199d6ff53..7e5bd8301 100644
--- a/docs/src/diffusions.md
+++ b/docs/src/diffusions.md
@@ -44,7 +44,7 @@ Or more compactly:
 | Time-fixed   | [`FixedDiffusion`](@ref)   | [`FixedMVDiffusion`](@ref)            |
 
 
-For more details on diffusions and calibration, check out this paper [[1]](@ref diffusionrefs).
+For more details on diffusions and calibration, check out this paper [bosch20capos](@cite).
 
 
 ## API
@@ -59,4 +59,9 @@ FixedMVDiffusion
 
 ## [References](@id diffusionrefs)
 
-[1] N. Bosch, P. Hennig, F. Tronarp: **Calibrated Adaptive Probabilistic ODE Solvers** (2021) ([link](http://proceedings.mlr.press/v130/bosch21a.html))
+```@bibliography
+Pages = []
+Canonical = false
+
+bosch20capos
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index 392a7762f..69fcfb000 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -37,41 +37,6 @@ For a quick introduction check out the "[Solving ODEs with Probabilistic Numeric
 - Compatible with DAEs in mass-matrix ODE form (see [Solving DAEs with Probabilistic Numerics](@ref))
 
 
-## [References](@id references)
-
-The main references _for this package_ include:
-
-- N. Bosch, P. Hennig, F. Tronarp:
-  [**Probabilistic Exponential Integrators**](https://arxiv.org/abs/2305.14978)
-  (2023)
-- N. Bosch, F. Tronarp, P. Hennig:
-  [**Pick-and-Mix Information Operators for Probabilistic ODE Solvers**](https://proceedings.mlr.press/v151/bosch22a.html)
-  (2022)
-- N. Krämer, N. Bosch, J. Schmidt, P. Hennig:
-  [**Probabilistic ODE Solutions in Millions of Dimensions**](https://proceedings.mlr.press/v162/kramer22b.html)
-  (2022)
-- N. Bosch, P. Hennig, F. Tronarp:
-  [**Calibrated Adaptive Probabilistic ODE Solvers**](http://proceedings.mlr.press/v130/bosch21a.html)
-  (2021)
-- F. Tronarp, S. Särkkä, and P. Hennig:
-  [**Bayesian ODE Solvers: The Maximum A Posteriori Estimate**](https://link.springer.com/article/10.1007/s11222-021-09993-7)
-  (2021)
-- N. Krämer, P. Hennig:
-  [**Stable Implementation of Probabilistic ODE Solvers**](https://arxiv.org/abs/2012.10106v1)
-  (2020)
-- H. Kersting, T. J. Sullivan, and P. Hennig:
-  [**Convergence Rates of Gaussian Ode Filters**](https://link.springer.com/article/10.1007/s11222-020-09972-4)
-  (2020)
-- F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig:
-  [**Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective**](https://link.springer.com/article/10.1007/s11222-019-09900-1)
-  (2019)
-- M. Schober, S. Särkkä, and P. Hennig:
-  [**A Probabilistic Model for the Numerical Solution of Initial Value Problems**](https://link.springer.com/article/10.1007/s11222-017-9798-7)
-  (2018)
-
-More references on ODE filters and on probabilistic numerics in general can be found on [probabilistic-numerics.org ](https://www.probabilistic-numerics.org/research/general/).
-
-
 ## Related packages
 
 - [probdiffeq](https://pnkraemer.github.io/probdiffeq/): Fast and feature-rich filtering-based probabilistic ODE solvers in JAX.
diff --git a/docs/src/initialization.md b/docs/src/initialization.md
index 3ed5f2bc7..7513f8e89 100644
--- a/docs/src/initialization.md
+++ b/docs/src/initialization.md
@@ -65,7 +65,13 @@ TaylorModeInit
 ClassicSolverInit
 ```
 
+
 ## [References](@id initrefs)
 
-[1] N. Krämer, P. Hennig: **Stable Implementation of Probabilistic ODE Solvers** (2020) ([link](https://arxiv.org/abs/2012.10106))
-[2] M. Schober, S. Särkkä, and P. Hennig: **A Probabilistic Model for the Numerical Solution of Initial Value Problems** (2018) ([link](https://link.springer.com/article/10.1007/s11222-017-9798-7))
+```@bibliography
+Pages = []
+Canonical = false
+
+kraemer20stableimplementation
+schober16probivp
+```
diff --git a/docs/src/references.md b/docs/src/references.md
new file mode 100644
index 000000000..50b6a16dc
--- /dev/null
+++ b/docs/src/references.md
@@ -0,0 +1,5 @@
+# [References](@id references)
+
+```@bibliography
+*
+```
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
new file mode 100644
index 000000000..136453adc
--- /dev/null
+++ b/docs/src/refs.bib
@@ -0,0 +1,172 @@
+@misc{bosch23expint,
+  title =        {Probabilistic Exponential Integrators},
+  author =       {Nathanael Bosch and Philipp Hennig and Filip Tronarp},
+  year =         2023,
+  eprint =       {2305.14978},
+  archivePrefix ={arXiv},
+  primaryClass = {math.NA}
+}
+
+@InProceedings{krämer21highdim,
+  title =        {Probabilistic {ODE} Solutions in Millions of Dimensions},
+  author =       {Krämer, Nicholas and Bosch, Nathanael and Schmidt, Jonathan
+                  and Hennig, Philipp},
+  booktitle =    {Proceedings of the 39th International Conference on Machine
+                  Learning},
+  pages =        {11634--11649},
+  year =         2022,
+  editor =       {Chaudhuri, Kamalika and Jegelka, Stefanie and Song, Le and
+                  Szepesvari, Csaba and Niu, Gang and Sabato, Sivan},
+  volume =       162,
+  series =       {Proceedings of Machine Learning Research},
+  month =        {17--23 Jul},
+  publisher =    {PMLR},
+  pdf =          {https://proceedings.mlr.press/v162/kramer22b/kramer22b.pdf},
+  url =          {https://proceedings.mlr.press/v162/kramer22b.html},
+}
+
+@InProceedings{bosch22pickandmix,
+  title =        {Pick-and-Mix Information Operators for Probabilistic {ODE}
+                  Solvers},
+  author =       {Bosch, Nathanael and Tronarp, Filip and Hennig, Philipp},
+  booktitle =    {Proceedings of The 25th International Conference on Artificial
+                  Intelligence and Statistics},
+  pages =        {10015--10027},
+  year =         2022,
+  editor =       {Camps-Valls, Gustau and Ruiz, Francisco J. R. and Valera,
+                  Isabel},
+  volume =       151,
+  series =       {Proceedings of Machine Learning Research},
+  month =        {28--30 Mar},
+  publisher =    {PMLR},
+  pdf =          {https://proceedings.mlr.press/v151/bosch22a/bosch22a.pdf},
+  url =          {https://proceedings.mlr.press/v151/bosch22a.html},
+}
+
+@InProceedings{bosch20capos,
+  title =        {Calibrated Adaptive Probabilistic {ODE} Solvers},
+  author =       {Bosch, Nathanael and Hennig, Philipp and Tronarp, Filip},
+  booktitle =    {Proceedings of The 24th International Conference on Artificial
+                  Intelligence and Statistics},
+  pages =        {3466--3474},
+  year =         2021,
+  editor =       {Banerjee, Arindam and Fukumizu, Kenji},
+  volume =       130,
+  series =       {Proceedings of Machine Learning Research},
+  month =        {13--15 Apr},
+  publisher =    {PMLR},
+  pdf =          {http://proceedings.mlr.press/v130/bosch21a/bosch21a.pdf},
+  url =          {http://proceedings.mlr.press/v130/bosch21a.html},
+}
+
+@article{tronarp20bayes,
+  author =       {Tronarp, Filip and Särkkä, Simo and Hennig, Philipp},
+  title =        {Bayesian {ODE} solvers: the maximum a posteriori estimate},
+  journal =      {Statistics and Computing},
+  year =         2021,
+  month =        {Mar},
+  day =          03,
+  volume =       31,
+  number =       3,
+  pages =        23,
+  issn =         {1573-1375},
+  doi =          {10.1007/s11222-021-09993-7},
+  url =          {https://doi.org/10.1007/s11222-021-09993-7}
+}
+
+@article{tronarp18probsol,
+  author =       {Filip Tronarp and Hans Kersting and Simo Särkkä and Philipp
+                  Hennig},
+  title =        {Probabilistic solutions to ordinary differential equations as
+                  nonlinear {B}ayesian filtering: a new perspective},
+  year =         2019,
+  volume =       29,
+  number =       6,
+  pages =        {1297-1315},
+  doi =          {10.1007/s11222-019-09900-1},
+  url =          {https://doi.org/10.1007/s11222-019-09900-1},
+  journal =      {Statistics and Computing},
+}
+
+@article{kraemer20stableimplementation,
+  author =       {Krämer, Nicholas and Hennig, Philipp},
+  title =        {Stable Implementation of Probabilistic {ODE} Solvers},
+  journal =      {CoRR},
+  year =         2020,
+  url =          {http://arxiv.org/abs/2012.10106v1},
+  archivePrefix ={arXiv},
+  eprint =       {2012.10106},
+  primaryClass = {stat.ML},
+}
+
+@article{kersting18convergencerates,
+  author =       {Kersting, Hans and Sullivan, T. J. and Hennig, Philipp},
+  title =        {Convergence rates of {G}aussian {ODE} filters},
+  journal =      {Statistics and Computing},
+  year =         2020,
+  month =        {Nov},
+  day =          01,
+  volume =       30,
+  number =       6,
+  pages =        {1791-1816},
+  issn =         {1573-1375},
+  doi =          {10.1007/s11222-020-09972-4},
+  url =          {https://doi.org/10.1007/s11222-020-09972-4}
+}
+
+@article{schober16probivp,
+  author =       "Schober, Michael and Särkkä, Simo and Hennig, Philipp",
+  title =        "A probabilistic model for the numerical solution of initial
+                  value problems",
+  journal =      "Statistics and Computing",
+  year =         2019,
+  month =        "Jan",
+  day =          01,
+  volume =       29,
+  number =       1,
+  pages =        "99--122",
+  issn =         "1573-1375",
+  doi =          "10.1007/s11222-017-9798-7",
+  url =          "https://doi.org/10.1007/s11222-017-9798-7"
+}
+
+@InProceedings{tronarp22fenrir,
+  title =        {Fenrir: Physics-Enhanced Regression for Initial Value
+                  Problems},
+  author =       {Tronarp, Filip and Bosch, Nathanael and Hennig, Philipp},
+  booktitle =    {Proceedings of the 39th International Conference on Machine
+                  Learning},
+  pages =        {21776--21794},
+  year =         2022,
+  editor =       {Chaudhuri, Kamalika and Jegelka, Stefanie and Song, Le and
+                  Szepesvari, Csaba and Niu, Gang and Sabato, Sivan},
+  volume =       162,
+  series =       {Proceedings of Machine Learning Research},
+  month =        {17--23 Jul},
+  publisher =    {PMLR},
+  pdf =          {https://proceedings.mlr.press/v162/tronarp22a/tronarp22a.pdf},
+  url =          {https://proceedings.mlr.press/v162/tronarp22a.html},
+}
+
+@inproceedings{kersting16active,
+  author =       {Kersting, Hans and Hennig, Philipp},
+  title =        {Active Uncertainty Calibration in Bayesian ODE Solvers},
+  year =         2016,
+  isbn =         9780996643115,
+  publisher =    {AUAI Press},
+  booktitle =    {Proceedings of the Thirty-Second Conference on Uncertainty in
+                  Artificial Intelligence},
+  pages =        {309–318},
+  numpages =     10,
+  series =       {UAI'16},
+  url =          {http://www.auai.org/uai2016/proceedings/papers/163.pdf},
+}
+
+@book{hennig22probnum,
+  place =        {Cambridge},
+  title =        {Probabilistic Numerics: Computation as Machine Learning},
+  DOI =          {10.1017/9781316681411},
+  publisher =    {Cambridge University Press},
+  author =       {Hennig, Philipp and Osborne, Michael A. and Kersting, Hans P.},
+  year =         2022,
+}
diff --git a/docs/src/solvers.md b/docs/src/solvers.md
index 62eb5844c..3dc5abe00 100644
--- a/docs/src/solvers.md
+++ b/docs/src/solvers.md
@@ -21,4 +21,10 @@ RosenbrockExpEK
 
 ## [References](@id solversrefs)
 
-[1] N. Bosch, F. Tronarp, P. Hennig: **Pick-and-Mix Information Operators for Probabilistic ODE Solvers** (2022) ([link](https://proceedings.mlr.press/v151/bosch22a.html))
+
+```@bibliography
+Pages = []
+Canonical = false
+
+bosch23expint
+```
diff --git a/docs/src/tutorials/dae.md b/docs/src/tutorials/dae.md
index 8152ef342..2e61e4cbf 100644
--- a/docs/src/tutorials/dae.md
+++ b/docs/src/tutorials/dae.md
@@ -136,4 +136,9 @@ So it seems that, even if the index-3 DAE could also be solved directly, index l
 
 ### References
 
-[1] N. Bosch, F. Tronarp, P. Hennig: **Pick-and-Mix Information Operators for Probabilistic ODE Solvers** (2022) ([link](https://proceedings.mlr.press/v151/bosch22a.html))
+```@bibliography
+Pages = []
+Canonical = false
+
+bosch22pickandmix
+```
diff --git a/docs/src/tutorials/dynamical_odes.md b/docs/src/tutorials/dynamical_odes.md
index 497dcbf07..b8f0655a9 100644
--- a/docs/src/tutorials/dynamical_odes.md
+++ b/docs/src/tutorials/dynamical_odes.md
@@ -85,7 +85,7 @@ plot(sol2, vars=(3, 4))
 Solving second-order ODEs is not just a matter of convenience - in fact, SciMLBase's `SecondOrderODEProblem` is neatly designed in such a way that all the classic solvers from OrdinaryDiffEq.jl can handle it by solving the corresponding first-order ODE.
 But, transforming the ODE to first order increases the dimensionality of the problem, and comes therefore at increased computational cost; this also motivates [classic specialized solvers for second-order ODEs](https://diffeq.sciml.ai/stable/solvers/dynamical_solve/).
 
-The probabilistic numerical solvers from ProbNumDiffEq.jl have the same internal state representation for first and second order ODEs; all that changes is the _measurement model_ [1].
+The probabilistic numerical solvers from ProbNumDiffEq.jl have the same internal state representation for first and second order ODEs; all that changes is the _measurement model_ [bosch22pickandmix](@citep).
 As a result, we can use the `EK1` both for first and second order ODEs, but it automatically specializes on the latter to provide a __2x performance boost__:
 
 ```
@@ -131,7 +131,7 @@ Let's fix this to get a physically more meaningful solution.
 
 ### Energy preservation with the `ManifoldUpdate` callback
 
-In the language of ODE filters, preserving energy over time amounts to just another measurement model [1].
+In the language of ODE filters, preserving energy over time amounts to just another measurement model [bosch22pickandmix](@citep).
 The most convenient way of updating on this additional zero measurement with ProbNumDiffEq.jl is with the `ManifoldUpdate` callback.
 
 !!! note
@@ -155,6 +155,12 @@ plot!(longsol_preserving.t, E.(longsol_preserving.u))
 
 Voilà! With the `ManifoldUpdate` callback we could preserve the energy over time and obtain a more truthful probabilistic numerical long-term simulation of the Hénon-Heiles model.
 
-#### References
 
-[1] N. Bosch, F. Tronarp, P. Hennig: **Pick-and-Mix Information Operators for Probabilistic ODE Solvers** (2022)
+### References
+
+```@bibliography
+Pages = []
+Canonical = false
+
+bosch22pickandmix
+```
diff --git a/docs/src/tutorials/exponential_integrators.md b/docs/src/tutorials/exponential_integrators.md
index 033e43c04..b23059aa9 100644
--- a/docs/src/tutorials/exponential_integrators.md
+++ b/docs/src/tutorials/exponential_integrators.md
@@ -8,7 +8,7 @@ Exponential integrators are a class of numerical methods for solving semi-linear
 ```
 where $L$ is a linear operator and $f$ is a nonlinear function.
 In a nutshell, exponential integrators solve the linear part of the ODE exactly, and only approximate the nonlinear part.
-[Probabilistic exponential integrators](https://arxiv.org/abs/2305.14978) are the probabilistic numerics approach to exponential integrators.
+[Probabilistic exponential integrators](@cite bosch23expint) [bosch23expint](@cite) are the probabilistic numerics approach to exponential integrators.
 
 ## Example
 
@@ -92,4 +92,9 @@ plot!(sol_expek1, denseplot=false, marker=:o, markersize=2, label="EK1 + IOUP")
 
 ### References
 
-[1] N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023) ([link](https://arxiv.org/abs/2305.14978))
+```@bibliography
+Pages = []
+Canonical = false
+
+bosch23expint
+```
diff --git a/docs/src/tutorials/fenrir.md b/docs/src/tutorials/fenrir.md
index 6db123e6b..d4f11feca 100644
--- a/docs/src/tutorials/fenrir.md
+++ b/docs/src/tutorials/fenrir.md
@@ -71,7 +71,7 @@ nll
 ```
 This is the negative marginal log-likelihood of the parameter `p_est`.
 You can use it as any other NLL: Optimize it to compute maximum-likelihood estimates or MAPs, or plug it into MCMC to sample from the posterior.
-In [our paper](https://proceedings.mlr.press/v162/tronarp22a.html) we compute MLEs by pairing Fenrir with [Optimization.jl](http://optimization.sciml.ai/stable/) and [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/).
+In our paper [tronarp22fenrir](@cite) we compute MLEs by pairing Fenrir with [Optimization.jl](http://optimization.sciml.ai/stable/) and [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/).
 Let's quickly explore how to do this next.
 
 
@@ -115,6 +115,11 @@ plot!(mle_sol, color=3, label=["MLE-parameter Solution" ""])
 Looks good!
 
 
-#### Reference
+### Reference
 
-[1] F. Tronarp, N. Bosch, P. Hennig: **Fenrir: Fenrir: Physics-Enhanced Regression for Initial Value Problems** [(2022)](https://proceedings.mlr.press/v162/tronarp22a.html)
+```@bibliography
+Pages = []
+Canonical = false
+
+tronarp22fenrir
+```
diff --git a/expdev.jl b/expdev.jl
new file mode 100644
index 000000000..bef5ad92e
--- /dev/null
+++ b/expdev.jl
@@ -0,0 +1,24 @@
+using ProbNumDiffEq, Plots, LinearAlgebra
+theme(:default; palette=["#4063D8", "#389826", "#9558B2", "#CB3C33"])
+
+f(du, u, p, t) = (@. du = -u + sin(u))
+u0 = [1.0]
+tspan = (0.0, 20.0)
+prob = ODEProblem(f, u0, tspan)
+
+ref = solve(prob, EK1(), abstol=1e-10, reltol=1e-10)
+
+STEPSIZE = 4
+DM = FixedDiffusion() # recommended for fixed steps
+
+sol_ek0 = solve(prob, EK0(prior=IOUP(3, -1), diffusionmodel=DM), adaptive=false,
+    dt=STEPSIZE)
+sol_ekl = solve(prob, EK0(prior=IOUP(3, -1), diffusionmodel=DM), adaptive=false,
+    dt=STEPSIZE)
+sol_ek1 = solve(prob, EK1(prior=IOUP(3, -1), diffusionmodel=DM), adaptive=false,
+    dt=STEPSIZE)
+
+plot(ref, color=:black, linestyle=:dash, label="Reference", ylims=(0.3, 1.05))
+plot!(sol_ek0, denseplot=false, marker=:o, markersize=2, label="EK0", color=1)
+plot!(sol_ekl, denseplot=false, marker=:o, markersize=2, label="EKL", color=2)
+plot!(sol_ek1, denseplot=false, marker=:o, markersize=2, label="EK1", color=3)
diff --git a/src/algorithms.jl b/src/algorithms.jl
index eb53e6251..efba620de 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -139,7 +139,7 @@ julia> solve(prob, ExpEK(L=-1))
 
 
 # Reference
-[1] N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023) ([link](https://arxiv.org/abs/2305.14978))
+* [bosch23expint](@cite) Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
 """
 ExpEK(; L, order=3, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)
 
@@ -172,7 +172,7 @@ julia> solve(prob, RosenbrockExpEK())
 ```
 
 # Reference
-[1] N. Bosch, P. Hennig, F. Tronarp: **Probabilistic Exponential Integrators** (2023) ([link](https://arxiv.org/abs/2305.14978))
+* [bosch23expint](@cite) Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
 """
 RosenbrockExpEK(; order=3, kwargs...) =
     EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)
diff --git a/src/diffusions.jl b/src/diffusions.jl
index 46a43350d..d3a600214 100644
--- a/src/diffusions.jl
+++ b/src/diffusions.jl
@@ -35,7 +35,8 @@ a diagonal matrix is estimated. This can be helpful to get more expressive poste
 covariances when using the [`EK0`](@ref), since the individual dimensions can be adjusted
 separately.
 
-See also [[1]](@ref diffusionrefs).
+# References
+* [bosch20capos](@cite) Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
 """
 struct DynamicMVDiffusion <: AbstractDynamicDiffusion end
 initial_diffusion(::DynamicMVDiffusion, d, q, Eltype) =
@@ -100,7 +101,8 @@ a diagonal matrix is estimated. This can be helpful to get more expressive poste
 covariances when using the [`EK0`](@ref), since the individual dimensions can be adjusted
 separately.
 
-See also [[1]](@ref diffusionrefs).
+# References
+* [bosch20capos](@cite) Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
 """
 Base.@kwdef struct FixedMVDiffusion{T} <: AbstractStaticDiffusion
     initial_diffusion::T = 1.0
@@ -149,7 +151,7 @@ Corresponds to
 where ``z, H, Q`` are taken from the passed integrator.
 
 For more background information
-- N. Bosch, P. Hennig, F. Tronarp: **Calibrated Adaptive Probabilistic ODE Solvers** (2021)
+* [bosch20capos](@cite) Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
 """
 function local_scalar_diffusion(cache)
     @unpack d, R, H, Qh, measurement, m_tmp, Smat = cache
@@ -177,7 +179,7 @@ where ``z, H, Q`` are taken from the passed integrator.
 **This should only be used with the EK0!**
 
 For more background information
-- N. Bosch, P. Hennig, F. Tronarp: **Calibrated Adaptive Probabilistic ODE Solvers** (2021)
+* [bosch20capos](@cite) Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
 """
 function local_diagonal_diffusion(cache)
     @unpack d, q, H, Qh, measurement, m_tmp, tmp = cache
diff --git a/src/initialization/common.jl b/src/initialization/common.jl
index b9180e14e..8223b296e 100644
--- a/src/initialization/common.jl
+++ b/src/initialization/common.jl
@@ -16,7 +16,7 @@ given problem (typically because the problem definition does not allow for eleme
  `Taylor`). If this happens, try [`ClassicSolverInit`](@ref).
 
 # References
-[[1]](@ref initrefs) N. Krämer, P. Hennig: **Stable Implementation of Probabilistic ODE Solvers** (2020)
+* [kraemer20stableimplementation](@cite) Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
 """
 struct TaylorModeInit <: InitializationScheme end
 
@@ -41,8 +41,8 @@ optionally the second derivative can also be set via automatic differentiation b
   automatic differentiation with ForwardDiff.jl.
 
 # References
-[[2]](@ref initrefs) M. Schober, S. Särkkä, and P. Hennig: **A Probabilistic Model for the Numerical Solution of Initial Value Problems** (2018)
-
+* [kraemer20stableimplementation](@cite) Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
+* [schober16probivp](@cite) Schober et al, "A probabilistic model for the numerical solution of initial value problems", Statistics and Computing (2019)
 """
 Base.@kwdef struct ClassicSolverInit{ALG} <: InitializationScheme
     alg::ALG = Tsit5()

From 7e471e84629bef880d56fbaa3767df8ab5b9dd4c Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Tue, 1 Aug 2023 09:24:35 +0000
Subject: [PATCH 9/9] Set version to 0.12.1

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 92e2366ba..6b3ec4bd3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ProbNumDiffEq"
 uuid = "bf3e78b0-7d74-48a5-b855-9609533b56a5"
 authors = ["Nathanael Bosch"]
-version = "0.12.0"
+version = "0.12.1"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"