diff --git a/src/DiffEqCallbacks.jl b/src/DiffEqCallbacks.jl
index 407fa71e..950bc26b 100644
--- a/src/DiffEqCallbacks.jl
+++ b/src/DiffEqCallbacks.jl
@@ -5,7 +5,6 @@ using DataStructures: DataStructures, BinaryMaxHeap, BinaryMinHeap
 using DiffEqBase: get_tstops, get_tstops_array, get_tstops_max
 using DifferentiationInterface: DifferentiationInterface, Constant
 using LinearAlgebra: LinearAlgebra, adjoint, axpy!, copyto!, diagind, mul!
-using Markdown: @doc_str
 using PrecompileTools: PrecompileTools
 using RecipesBase: @recipe
 using RecursiveArrayTools: RecursiveArrayTools, DiffEqArray, copyat_or_push!
diff --git a/src/autoabstol.jl b/src/autoabstol.jl
index 5ff19037..9280b152 100644
--- a/src/autoabstol.jl
+++ b/src/autoabstol.jl
@@ -42,6 +42,23 @@ is set to the maximum value that the state has thus far reached times the relati
 
 If this callback is used in isolation, `save=true` is required for normal saving behavior.
 Otherwise, `save=false` should be set to ensure extra saves do not occur.
+
+## Returns
+
+A `DiscreteCallback` that updates the integrator absolute tolerance from the largest
+state magnitude observed so far.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+f(u, p, t) = 0.5u
+prob = ODEProblem(f, 1.0, (0.0, 2.0))
+cb = AutoAbstol(; init_curmax = 1.0e-8)
+
+sol = solve(prob, Tsit5(); callback = cb, reltol = 1.0e-6)
+```
 """
 function AutoAbstol(save = true; init_curmax = 0.0)
     affect! = AutoAbstolAffect(abs.(init_curmax))
diff --git a/src/domain.jl b/src/domain.jl
index 33fb19aa..42a52244 100644
--- a/src/domain.jl
+++ b/src/domain.jl
@@ -280,6 +280,20 @@ inside the domain. Thus, a `PositiveDomain` callback should generally be preferr
 Shampine, Lawrence F., Skip Thompson, Jacek Kierzenka and G. D. Byrne.
 Non-negative solutions of ODEs. Applied Mathematics and Computation 170
 (2005): 556-569.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+function nonnegative_residual(resid, u, p, t)
+    @. resid = max(-u, 0)
+end
+
+prob = ODEProblem((du, u, p, t) -> (du .= -u), [1.0, 2.0], (0.0, 2.0))
+cb = GeneralDomain(nonnegative_residual, [1.0, 2.0]; abstol = 1.0e-8)
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function GeneralDomain(
         g, u = nothing; save = true, abstol = nothing, scalefactor = nothing,
@@ -310,15 +324,13 @@ function GeneralDomain(
     return CallbackSet(manifold_projection, domain_cb)
 end
 
-@doc doc"""
-```julia
-PositiveDomain(u = nothing; save = true, abstol = nothing, scalefactor = nothing)
-```
+"""
+    PositiveDomain(u = nothing; save = true, abstol = nothing, scalefactor = nothing)
 
 Especially in biology and other natural sciences, a desired property of
 dynamical systems is the positive invariance of the positive cone, i.e.
 non-negativity of variables at time ``t_0`` ensures their non-negativity at times
-``t \geq t_0`` for which the solution is defined. However, even if a system
+``t \\geq t_0`` for which the solution is defined. However, even if a system
 satisfies this property mathematically it can be difficult for ODE solvers to
 ensure it numerically, as these [MATLAB examples](https://www.mathworks.com/help/matlab/math/nonnegative-ode-solution.html)
 show.
@@ -346,7 +358,7 @@ depends on how accurately extrapolations approximate next time steps.
 Please note, that the system should be defined also outside the positive domain,
 since even with these approaches, negative variables might occur during the
 calculations. Moreover, one should follow Shampine's et al. advice and set the
-derivative ``x'_i`` of a negative component ``x_i`` to ``\max \{0, f_i(x, t)\}``,
+derivative ``x'_i`` of a negative component ``x_i`` to ``\\max \\{0, f_i(x, t)\\}``,
 where ``t`` denotes the current time point with state vector ``x`` and ``f_i``
 is the ``i``-th component of function ``f`` in an ODE system ``x' = f(x, t)``.
 
@@ -370,6 +382,18 @@ is the ``i``-th component of function ``f`` in an ODE system ``x' = f(x, t)``.
 Shampine, Lawrence F., Skip Thompson, Jacek Kierzenka and G. D. Byrne.
 Non-negative solutions of ODEs. Applied Mathematics and Computation 170
 (2005): 556-569.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+f(u, p, t) = -u
+prob = ODEProblem(f, 1.0, (0.0, 2.0))
+cb = PositiveDomain()
+
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function PositiveDomain(u = nothing; save = true, abstol = nothing, scalefactor = nothing)
     if u isa Nothing
diff --git a/src/function_caller.jl b/src/function_caller.jl
index 32efd7cf..e8b20f49 100644
--- a/src/function_caller.jl
+++ b/src/function_caller.jl
@@ -60,23 +60,43 @@ function functioncalling_initialize(cb, u, t, integrator)
 end
 
 """
-```julia
-FunctionCallingCallback(func;
-    funcat = Vector{Float64}(),
-    func_everystep = isempty(funcat),
-    func_start = true,
-    tdir = 1)
-```
+    FunctionCallingCallback(func;
+        funcat = Vector{Float64}(),
+        func_everystep = isempty(funcat),
+        func_start = true,
+        tdir = 1)
 
 The function calling callback lets you define a function `func(u,t,integrator)`
-which gets called at the time points of interest. The constructor is:
+which gets called at the time points of interest.
+
+## Arguments
 
   - `func(u, t, integrator)` is the function to be called.
+
+## Keyword Arguments
+
   - `funcat` values or interval that the function is sure to be evaluated at.
   - `func_everystep` whether to call the function after each integrator step.
   - `func_start` whether the function is called at the initial condition.
   - `tdir` should be `sign(tspan[end]-tspan[1])`. It defaults to `1` and should
     be adapted if `tspan[1] > tspan[end]`.
+
+## Returns
+
+A `DiscreteCallback` that calls `func` without modifying the integrator state.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+seen = Float64[]
+func = (u, t, integrator) -> push!(seen, t)
+cb = FunctionCallingCallback(func; funcat = 0.0:0.25:1.0)
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function FunctionCallingCallback(
         func;
diff --git a/src/integrating.jl b/src/integrating.jl
index 770dc2c0..289cfb0e 100644
--- a/src/integrating.jl
+++ b/src/integrating.jl
@@ -235,6 +235,23 @@ returns Integral(integrand_func(u(t),t)dt over the problem tspan.
 The outputted values are saved into `integrand_values`. The values are found
 via `integrand_values.integrand`.
 
+## Returns
+
+A `DiscreteCallback` that saves one quadrature estimate per accepted step.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+values = IntegrandValues(Float64, Float64)
+cb = IntegratingCallback((u, t, integrator) -> u^2, values, 0.0)
+
+sol = solve(prob, Tsit5(); callback = cb)
+integrals_by_step = values.integrand
+```
+
 !!! note
 
     This method is currently limited to ODE solvers of order 10 or lower. Open an issue if other
diff --git a/src/integrating_GK_affect.jl b/src/integrating_GK_affect.jl
index dc2c1ad8..6fb811cb 100644
--- a/src/integrating_GK_affect.jl
+++ b/src/integrating_GK_affect.jl
@@ -269,6 +269,24 @@ returns Integral(integrand_func(u(t),t)dt) over the problem tspan.
 The outputted values are saved into `integrand_values`. The values are found
 via `integrand_values.integrand`.
 
+## Returns
+
+A `DiscreteCallback` that saves one Gauss-Kronrod quadrature estimate per
+accepted step.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+values = IntegrandValues(Float64, Float64)
+cb = IntegratingGKCallback((u, t, integrator) -> u^2, values, 0.0)
+
+sol = solve(prob, Tsit5(); callback = cb)
+integrals_by_step = values.integrand
+```
+
 !!! note
 
     Method has automatic error control (h-adaptive quadrature).
diff --git a/src/integrating_GK_sum.jl b/src/integrating_GK_sum.jl
index 5ebbcdcd..17ca6d80 100644
--- a/src/integrating_GK_sum.jl
+++ b/src/integrating_GK_sum.jl
@@ -85,9 +85,11 @@ end
 
 """
 ```julia
-IntegratingCallback(integrand_func,
-    integrand_values::IntegrandValues,
-    cache = nothing)
+IntegratingGKSumCallback(integrand_func,
+    integrand_values::IntegrandValuesSum,
+    integrand_prototype,
+    tol = 1.0e-7;
+    integrand_inplace = nothing)
 ```
 
 Lets one define a function `integrand_func(u, t, integrator)` which
@@ -119,6 +121,24 @@ returns Integral(integrand_func(u(t),t)dt over the problem tspan.
 The outputted values are saved into `integrand_values`. The values are found
 via `integrand_values.integrand`.
 
+## Returns
+
+A `DiscreteCallback` that accumulates Gauss-Kronrod quadrature estimates into
+`integrand_values.integrand`.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+values = IntegrandValuesSum(0.0)
+cb = IntegratingGKSumCallback((u, t, integrator) -> u^2, values, 0.0)
+
+sol = solve(prob, Tsit5(); callback = cb)
+total = values.integrand
+```
+
 !!! note
 
     This method uses Gauss-Kronrod quadrature rule to allow for error control.
diff --git a/src/integrating_sum.jl b/src/integrating_sum.jl
index 02d8bef2..fe5098b2 100644
--- a/src/integrating_sum.jl
+++ b/src/integrating_sum.jl
@@ -119,9 +119,10 @@ end
 
 """
 ```julia
-IntegratingCallback(integrand_func,
-    integrand_values::IntegrandValues,
-    cache = nothing)
+IntegratingSumCallback(integrand_func,
+    integrand_values::IntegrandValuesSum,
+    integrand_prototype;
+    integrand_inplace = nothing)
 ```
 
 Lets one define a function `integrand_func(u, t, integrator)` which
@@ -153,6 +154,24 @@ returns Integral(integrand_func(u(t),t)dt over the problem tspan.
 The outputted values are saved into `integrand_values`. The values are found
 via `integrand_values.integrand`.
 
+## Returns
+
+A `DiscreteCallback` that accumulates the quadrature estimate into
+`integrand_values.integrand`.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+values = IntegrandValuesSum(0.0)
+cb = IntegratingSumCallback((u, t, integrator) -> u^2, values, 0.0)
+
+sol = solve(prob, Tsit5(); callback = cb)
+total = values.integrand
+```
+
 !!! note
 
     This method is currently limited to ODE solvers of order 10 or lower. Open an issue if other
diff --git a/src/iterative_and_periodic.jl b/src/iterative_and_periodic.jl
index 423382d5..6c2dbd05 100644
--- a/src/iterative_and_periodic.jl
+++ b/src/iterative_and_periodic.jl
@@ -1,8 +1,6 @@
 """
-```julia
-IterativeCallback(time_choice, user_affect!, tType = Float64;
-    initial_affect = false, kwargs...)
-```
+    IterativeCallback(time_choice, user_affect!, tType = Float64;
+        initial_affect = false, kwargs...)
 
 A callback to be used to iteratively apply some affect. For example, if given the first
 effect at `t₁`, you can define `t₂` to apply the next effect.
@@ -16,6 +14,28 @@ effect at `t₁`, you can define `t₂` to apply the next effect.
 ## Keyword Arguments
 
   - `initial_affect` is whether to apply the affect at `t=0` which defaults to `false`
+  - `initialize` is the callback initialization function.
+  - `kwargs` are keyword arguments accepted by `DiscreteCallback`.
+
+## Returns
+
+A `DiscreteCallback` that repeatedly schedules the next time returned by
+`time_choice`.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+count = Ref(0)
+hits = Float64[]
+time_choice = integrator -> (count[] += 1; count[] <= 3 ? integrator.t + 0.1 : nothing)
+affect! = integrator -> push!(hits, integrator.t)
+cb = IterativeCallback(time_choice, affect!)
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function IterativeCallback(
         time_choice, user_affect!, tType = Float64;
@@ -120,11 +140,8 @@ function add_next_tstop!(integrator, S::PeriodicCallbackAffect)
 end
 
 """
-```julia
-PeriodicCallback(f, Δt::Number; phase = 0, initial_affect = false,
-    final_affect = false,
-    kwargs...)
-```
+    PeriodicCallback(f, Δt::Number; phase = 0, initial_affect = false,
+        final_affect = false, kwargs...)
 
 `PeriodicCallback` can be used when a function should be called periodically in terms of
 integration time (as opposed to wall time), i.e. at `t = tspan[1]`, `t = tspan[1] + Δt`,
@@ -149,6 +166,24 @@ discrete-time controller for a continuous-time system, running at a fixed rate.
   - `initial_affect` is whether to apply the affect at the initial time, which defaults to `false`
   - `final_affect` is whether to apply the affect at the final time, which defaults to `false`
   - `kwargs` are keyword arguments accepted by the `DiscreteCallback` constructor.
+
+## Returns
+
+A `DiscreteCallback` that schedules an `affect!` call every `Δt` units of
+integration time.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+samples = Float64[]
+affect! = integrator -> push!(samples, integrator.u)
+cb = PeriodicCallback(affect!, 0.1; initial_affect = true)
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function PeriodicCallback(
         f, Δt::Number;
diff --git a/src/manifold.jl b/src/manifold.jl
index 0382703d..a9dbd17d 100644
--- a/src/manifold.jl
+++ b/src/manifold.jl
@@ -66,6 +66,26 @@ order of the integrator even with the ManifoldProjection.
 [1] Ernst Hairer, Christian Lubich, Gerhard Wanner. Geometric Numerical Integration:
 Structure-Preserving Algorithms for Ordinary Differential Equations. Berlin ;
 New York :Springer, 2002.
+
+## Examples
+
+```julia
+using ADTypes, DiffEqCallbacks, OrdinaryDiffEq
+
+function unit_circle(resid, u, p, t)
+    resid[1] = sum(abs2, u) - 1
+end
+
+function rotation!(du, u, p, t)
+    du[1] = -u[2]
+    du[2] = u[1]
+end
+
+prob = ODEProblem(rotation!, [1.0, 0.0], (0.0, 10.0))
+cb = ManifoldProjection(unit_circle; residual_prototype = [0.0], autodiff = AutoFiniteDiff())
+
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 @concrete mutable struct ManifoldProjection
     manifold
diff --git a/src/preset_time.jl b/src/preset_time.jl
index 9bc8fbb1..17aaf971 100644
--- a/src/preset_time.jl
+++ b/src/preset_time.jl
@@ -32,12 +32,10 @@ function (f::PresetTimeFunction)(c, u, t, integrator)
 end
 
 """
-```julia
-PresetTimeCallback(tstops, user_affect!;
-    initialize = INITIALIZE_DEFAULT,
-    filter_tstops = true,
-    kwargs...)
-```
+    PresetTimeCallback(tstops, user_affect!;
+        initialize = INITIALIZE_DEFAULT,
+        filter_tstops = true,
+        kwargs...)
 
 A callback that adds callback `affect!` calls at preset times. No playing around with
 `tstops` or anything is required: this callback adds the triggers for you to make it
@@ -52,6 +50,26 @@ automatic.
 
   - `filter_tstops`: Whether to filter out tstops beyond the end of the integration timespan.
     Defaults to true. If false, then tstops can extend the interval of integration.
+  - `initialize`: callback initialization function.
+  - `kwargs`: keyword arguments forwarded to `DiscreteCallback`.
+
+## Returns
+
+A `DiscreteCallback` that schedules the requested `tstops` and calls `user_affect!`
+at those times.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+hits = Float64[]
+affect! = integrator -> push!(hits, integrator.t)
+cb = PresetTimeCallback([0.25, 0.5, 0.75], affect!)
+
+prob = ODEProblem((u, p, t) -> -u, 1.0, (0.0, 1.0))
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function PresetTimeCallback(
         tstops, user_affect!;
diff --git a/src/probints.jl b/src/probints.jl
index 32fc7ff0..630c406a 100644
--- a/src/probints.jl
+++ b/src/probints.jl
@@ -41,6 +41,17 @@ the timesteps and the order of the algorithm.
 
 Conrad P., Girolami M., Särkkä S., Stuart A., Zygalakis. K, Probability
 Measures for Numerical Solutions of Differential Equations, arXiv:1506.04592
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+prob = ODEProblem((u, p, t) -> -u, [1.0], (0.0, 1.0))
+cb = ProbIntsUncertainty(0.01, 5)
+
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function ProbIntsUncertainty(σ, order, save = true)
     affect! = ProbIntsCache(σ, order)
@@ -79,6 +90,17 @@ every step.
 
 Conrad P., Girolami M., Särkkä S., Stuart A., Zygalakis. K, Probability
 Measures for Numerical Solutions of Differential Equations, arXiv:1506.04592
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+prob = ODEProblem((u, p, t) -> -u, [1.0], (0.0, 1.0))
+cb = AdaptiveProbIntsUncertainty(5)
+
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function AdaptiveProbIntsUncertainty(order, save = true)
     affect! = AdaptiveProbIntsCache(order)
diff --git a/src/stepsizelimiters.jl b/src/stepsizelimiters.jl
index 542c3dc5..7f1b76bc 100644
--- a/src/stepsizelimiters.jl
+++ b/src/stepsizelimiters.jl
@@ -26,14 +26,13 @@ function StepsizeLimiter_initialize(cb, u, t, integrator)
     return cb.affect!(integrator)
 end
 
-@doc doc"""
-```julia
-StepsizeLimiter(dtFE;safety_factor=9//10,max_step=false,cached_dtcache=0.0)
-```
+"""
+    StepsizeLimiter(dtFE; safety_factor = 9 // 10, max_step = false,
+        cached_dtcache = 0.0)
 
 In many cases, there is a known maximal stepsize for which the computation is
 stable and produces correct results. For example, in hyperbolic PDEs one normally
-needs to ensure that the stepsize stays below some ``\Delta t_{FE}`` determined
+needs to ensure that the stepsize stays below some ``\\Delta t_{FE}`` determined
 by the CFL condition. For nonlinear hyperbolic PDEs this limit can be a function
 `dtFE(u,p,t)` which changes throughout the computation. The stepsize limiter lets
 you pass a function which will adaptively limit the stepsizes to match these
@@ -51,6 +50,23 @@ constraints.
   solver is set to `adaptive=false`.
 - `cached_dtcache` should be set to match the type for time when not using
   Float64 values.
+
+## Returns
+
+A `DiscreteCallback` that updates `integrator.opts.dtmax` before every step.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+f(u, p, t) = -u
+prob = ODEProblem(f, 1.0, (0.0, 1.0))
+dtFE(u, p, t) = 0.05
+cb = StepsizeLimiter(dtFE; safety_factor = 0.8)
+
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function StepsizeLimiter(
         dtFE; safety_factor = 9 // 10, max_step = false,
diff --git a/src/terminatesteadystate.jl b/src/terminatesteadystate.jl
index 91cca09c..8ea2e3ad 100644
--- a/src/terminatesteadystate.jl
+++ b/src/terminatesteadystate.jl
@@ -81,6 +81,22 @@ the [Steady State Solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/stead
 
   - `min_t` specifies an optional minimum `t` before the steady state calculations are allowed
     to terminate.
+
+## Returns
+
+A `DiscreteCallback` that terminates the integrator when `test` returns `true`.
+
+## Examples
+
+```julia
+using DiffEqCallbacks, OrdinaryDiffEq
+
+f(u, p, t) = 1 - u
+prob = ODEProblem(f, 0.0, (0.0, 100.0))
+cb = TerminateSteadyState(1.0e-8, 1.0e-8)
+
+sol = solve(prob, Tsit5(); callback = cb)
+```
 """
 function TerminateSteadyState(
         abstol = 1.0e-8, reltol = 1.0e-6, test::T = allDerivPass;
diff --git a/test/qa/qa.jl b/test/qa/qa.jl
index 82ff1c5d..74b96d4d 100644
--- a/test/qa/qa.jl
+++ b/test/qa/qa.jl
@@ -23,3 +23,40 @@ run_qa(
         ),
     ),
 )
+
+@testset "Public API documentation coverage" begin
+    public_names = Set(names(DiffEqCallbacks; all = false))
+    delete!(public_names, :DiffEqCallbacks)
+
+    doc_pages = joinpath(pkgdir(DiffEqCallbacks), "docs", "src")
+    rendered_names = Set{Symbol}()
+    for page in filter(endswith(".md"), readdir(doc_pages; join = true))
+        in_docs_block = false
+        for line in eachline(page)
+            stripped = strip(line)
+            if stripped == "```@docs"
+                in_docs_block = true
+                continue
+            elseif in_docs_block && startswith(stripped, "```")
+                in_docs_block = false
+                continue
+            end
+            if in_docs_block && !isempty(stripped) && !startswith(stripped, "#")
+                name = last(split(stripped, '.'))
+                push!(rendered_names, Symbol(name))
+            end
+        end
+    end
+
+    @testset "docstrings exist" begin
+        missing_docstrings = sort!(collect(filter(name -> !Docs.hasdoc(DiffEqCallbacks, name), public_names)))
+        @test isempty(missing_docstrings)
+    end
+
+    @testset "rendered docs entries exist" begin
+        missing_rendered_entries = sort!(collect(setdiff(public_names, rendered_names)))
+        stale_rendered_entries = sort!(collect(setdiff(rendered_names, public_names)))
+        @test isempty(missing_rendered_entries)
+        @test isempty(stale_rendered_entries)
+    end
+end