diff --git a/docs/src/examples/05f-cyclefree.jl b/docs/src/examples/05f-cyclefree.jl
index 78e7f922..68c66570 100644
--- a/docs/src/examples/05f-cyclefree.jl
+++ b/docs/src/examples/05f-cyclefree.jl
@@ -99,7 +99,7 @@ warmup = vcat(
             cs,
             cs.fluxes,
             optimizer = HiGHS.Optimizer,
-            output = (_, om) -> JuMP.value.(om[:x]),
+            output = (_, om) -> variable_vector(om),
             output_type = Vector{Float64},
         ) for v in vs
     )...,
@@ -114,7 +114,7 @@ sample = sample_constraints(
     seed = UInt(1234),
     output = cs.fluxes,
     n_chains = 10,
-    collect_iterations = collect(100:105),
+    collect_iterations = collect(10:15),
 )
 
 # The results can be observed (and usually plotted) from the sample vectors,
diff --git a/src/solver.jl b/src/solver.jl
index 2cc21d6f..38e58fdb 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -157,6 +157,16 @@ is_solved(opt_model::J.Model) =
 
 export is_solved
 
+"""
+$(TYPEDSIGNATURES)
+
+Retrieve the variable vector from a JuMP model created by
+[`optimization_model`](@ref).
+"""
+variable_vector(opt_model::J.Model) = J.value.(opt_model[:x])
+
+export variable_vector
+
 """
     Minimal