diff --git a/docs/literate/src/files/custom_semidiscretization.jl b/docs/literate/src/files/custom_semidiscretization.jl
index 609c76c0434..66f81d482ae 100644
--- a/docs/literate/src/files/custom_semidiscretization.jl
+++ b/docs/literate/src/files/custom_semidiscretization.jl
@@ -24,11 +24,125 @@
 # passed to [`semidiscretize`](@ref) as parameter.
 # For a [`SemidiscretizationHyperbolic`](@ref), the `ODEProblem` wraps
 # `Trixi.rhs!` as ODE RHS.
-#
+# For a [`SemidiscretizationHyperbolicParabolic`](@ref),  Trixi.jl
+# uses a `SplitODEProblem` combining `Trixi.rhs_parabolic!` for the
+# (potentially) stiff part and `Trixi.rhs!` for the other part.
+
+
+# ## Standard Trixi.jl setup
+
+# In this tutorial, we will consider the linear advection equation
+# with source term
+# ```math
+# \partial_t u(t,x) + \partial_x u(t,x) = -\exp(-t) sin\bigl(\pi (x - t) \bigr)
+# ```
+# with periodic boundary conditions in the domain `[-1, 1]` as a
+# model problem.
+# The initial condition is
+# ```math
+# u(0,x) = \sin(\pi x).
+# ```
+# The source term results in some damping and the analytical solution
+# ```math
+# u(t,x) = \exp(-t) \sin\bigl(\pi (x - t) \bigr).
+# ```
+# First, we discretize this equation using the standard functionality
+# of Trixi.jl.
+
+using Trixi, OrdinaryDiffEq, Plots
+
+# The linear scalar advection equation is already implemented in
+# Trixi.jl as [`LinearScalarAdvectionEquation1D`](@ref). We construct
+# it with an advection velocity `1.0`.
+
+equations = LinearScalarAdvectionEquation1D(1.0)
+
+# Next, we use a standard [`DGSEM`](@ref) solver.
+
+solver = DGSEM(polydeg = 3)
+
+# We create a simple [`TreeMesh`](@ref) in 1D.
+
+coordinates_min = (-1.0,)
+coordinates_max = (+1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max;
+                initial_refinement_level = 4,
+                n_cells_max = 10^4,
+                periodicity = true)
+
+# We wrap everything in in a semidiscretization and pass the source
+# terms as a standard Julia function. Please note that Trixi.jl uses
+# `SVector`s from
+# [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)
+# to store the conserved variables `u`. Thus, the return value of the
+# source terms must be wrapped in an `SVector` - even if we consider
+# just a scalar problem.
+
+function initial_condition(x, t, equations)
+    return SVector(exp(-t) * sinpi(x[1] - t))
+end
+
+function source_terms_standard(u, x, t, equations)
+    return -initial_condition(x, t, equations)
+end
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
+                                    solver;
+                                    source_terms = source_terms_standard)
+
+# Now, we can create the `ODEProblem`, solve the resulting ODE
+# using a time integration method from
+# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl),
+# and visualize the numerical solution at the final time using
+# [Plots.jl](https://github.com/JuliaPlots/Plots.jl).
+
+tspan = (0.0, 3.0)
+ode = semidiscretize(semi, tspan)
+
+sol = solve(ode, RDPK3SpFSAL49(); ode_default_options()...)
+
+plot(sol; label = "numerical sol.")
+
+# We can also plot the analytical solution for comparison.
+# Since Trixi.jl uses `SVector`s for the variables, we take their `first`
+# (and only) component to get the scalar value for manual plotting.
+
+let
+   x = range(-1.0, 1.0; length = 200)
+   plot!(x, first.(initial_condition.(x, sol.t[end], equations)),
+         label = "analytical sol.", linestyle = :dash, legend = :topleft)
+end
+
+# We can also add the initial condition to the plot.
+
+plot!(sol.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)
+
+# You can of course also use some
+# [callbacks](https://trixi-framework.github.io/Trixi.jl/stable/callbacks/)
+# provided by Trixi.jl as usual.
+
+summary_callback = SummaryCallback()
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi; interval = analysis_interval)
+alive_callback = AliveCallback(; analysis_interval)
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+sol = solve(ode, RDPK3SpFSAL49();
+            ode_default_options()..., callback = callbacks)
+summary_callback()
+
+
+# ## Using a custom ODE right-hand side function
+
+# TODO
 
 
 # ## Setting up a custom semidiscretization
 
+# TODO
+
 # Required
 # - `Trixi.rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t)`
 # - `Trixi.mesh_equations_solver_cache(semi::SemidiscretizationEulerGravity)`
diff --git a/docs/src/callbacks.md b/docs/src/callbacks.md
index 1d3e5e34b51..7e0ef127856 100644
--- a/docs/src/callbacks.md
+++ b/docs/src/callbacks.md
@@ -106,7 +106,7 @@ will yield the following plot:
   the automated performance measurements, including an output of the recorded timers after a simulation.
 * The [`VisualizationCallback`](@ref) can be used for in-situ visualization. See
   [Visualizing results during a simulation](@ref).
-* The [`TrivialCallback`](@ref) does nothing and can be used to to easily disable some callbacks
+* The [`TrivialCallback`](@ref) does nothing and can be used to easily disable some callbacks
   via [`trixi_include`](@ref).
 
 ### Equation-specific callbacks