diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
index 3af6956ac60..ec8ae43821f 100644
--- a/.github/workflows/Documenter.yml
+++ b/.github/workflows/Documenter.yml
@@ -3,7 +3,7 @@ name: Documentation
 on:
   push:
     branches:
-      - main
+      - 'main'
     tags: '*'
   pull_request:
 
@@ -21,4 +21,4 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
-        run: julia --project=docs/ docs/make.jl
+        run: julia --project=docs --color=yes docs/make.jl
diff --git a/docs/make.jl b/docs/make.jl
index 3595a87b8c6..02d01cd3c94 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -55,7 +55,8 @@ makedocs(
         "Authors" => "authors.md",
         "Contributing" => "contributing.md",
         "License" => "license.md"
-    ]
+    ],
+    strict = true # to make the GitHub action fail when doctests fail, see https://github.com/neuropsychology/Psycho.jl/issues/34
 )
 
 deploydocs(
diff --git a/docs/src/differentiable_programming.md b/docs/src/differentiable_programming.md
index 2f51cc577c1..89c0e839af5 100644
--- a/docs/src/differentiable_programming.md
+++ b/docs/src/differentiable_programming.md
@@ -78,11 +78,11 @@ julia> λ = eigvals(J);
 
 julia> scatter(real.(λ), imag.(λ));
 
-julia> round(maximum(real, λ) / maximum(abs, λ), sigdigits=2)
-6.7e-10
+julia> 3.0e-10 < maximum(real, λ) / maximum(abs, λ) < 8.0e-10
+true
 
-julia> round(maximum(real, λ), sigdigits=2)
-2.1e-7
+julia> 1.0e-7 < maximum(real, λ) < 5.0e-7
+true
 ```
 
 Interestingly, if we add dissipation by switching to the `flux_lax_friedrichs` at the interfaces,
@@ -133,16 +133,16 @@ julia> λ = eigvals(J);
 
 julia> scatter(real.(λ), imag.(λ));
 
-julia> round(maximum(real, λ) / maximum(abs, λ), sigdigits=2)
-3.2e-16
+julia> 1.0e-16 < maximum(real, λ) / maximum(abs, λ) < 6.0e-16
+true
 
-julia> round(maximum(real, λ), sigdigits=2)
-3.2e-12
+julia> 1.0e-12 < maximum(real, λ) < 6.0e-12
+true
 
 julia> λ, V = eigen(J);
 
-julia> round(cond(V), sigdigits=2)
-250.0
+julia> 200 < cond(V) < 300
+true
 ```
 
 If we add dissipation, the maximal real part is still approximately zero.
@@ -158,16 +158,18 @@ julia> λ = eigvals(J);
 
 julia> scatter!(real.(λ), imag.(λ));
 
-julia> λ = eigvals(J); round(maximum(real, λ) / maximum(abs, λ), sigdigits=2)
-5.3e-18
+julia> λ = eigvals(J);
 
-julia> round(maximum(real, λ), sigdigits=2)
-7.7e-14
+julia> 1.0e-18 < maximum(real, λ) / maximum(abs, λ) < 1.0e-16
+true
+
+julia> 5.0e-14 < maximum(real, λ) < 5.0e-13
+true
 
 julia> λ, V = eigen(J);
 
-julia> round(cond(V), sigdigits=2)
-93000.0
+julia> 90_000 < cond(V) < 100_000
+true
 ```
 Note that the condition number of the eigenvector matrix increases but is still smaller than for the
 example in 2D.