diff --git a/README.md b/README.md
index 789b3a8..bd26e4d 100644
--- a/README.md
+++ b/README.md
@@ -88,3 +88,18 @@ Currently, diagnostic functions for checking symmetry, symplecticity and the so-
 * `satisfies_simplifying_assumption_d(tab, σ=tab.s)`
 
 This list is expected to grow in the near future.
+
+
+## References
+
+If you use RungeKutta.jl in your work, please consider citing it by
+
+```
+@misc{Kraus:2020:RungeKutta,
+  title={RungeKutta.jl: Runge-Kutta Methods in Julia},
+  author={Kraus, Michael},
+  year={2020},
+  howpublished={\url{https://github.com/JuliaGNI/RungeKutta.jl}},
+  doi={10.5281/zenodo.4294923}
+}
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index 740eec7..7884e8e 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -64,3 +64,18 @@ In addition there exist functions to compute Gauss, Lobatto and Radau tableaus w
 
 All constructors take an optional type argument `T`, as in `TableauExplicitMidpoint(T)` or `TableauGauss(T,s)`. The default type is `Float64`, but it can be set to other number types if needed, e.g., to `Float32` for single precision or to the `Dec128` type from [DecFP.jl](https://github.com/JuliaMath/DecFP.jl) for quadruple precision.
 Internally, all tableaus are computed using `BigFloat`, providing high-accuracy coefficients as they are required for simulations in quadruple or higher precision. The internal precision can be set via `setprecision(40)`, cf. the [Julia Manual](https://docs.julialang.org/en/v1/) on [Arbitrary Precision Arithmetic](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic).
+
+
+## References
+
+If you use RungeKutta.jl in your work, please consider citing it by
+
+```
+@misc{Kraus:2020:RungeKutta,
+  title={RungeKutta.jl: Runge-Kutta Methods in Julia},
+  author={Kraus, Michael},
+  year={2020},
+  howpublished={\url{https://github.com/JuliaGNI/RungeKutta.jl}},
+  doi={10.5281/zenodo.4294923}
+}
+```