diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 23dd49f9..1449e203 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -53,6 +53,10 @@ jobs:
       - uses: julia-actions/setup-julia@v1
         with:
           version: '1'
+      - uses: julia-actions/cache@v1
+        with:
+          cache-registries: "true"
+          cache-compiled: "true"
       - name: Configure doc environment
         run: |
           julia --project=docs/ -e '
diff --git a/.gitignore b/.gitignore
index 3797f8d8..d2765d9b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@
 /docs/build/
 /test/Manifest.toml
 .vscode
+/docs/src/generated/*.md
diff --git a/docs/Project.toml b/docs/Project.toml
index b1dffb06..bb9a92a8 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,5 +4,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 HMatrices = "8646bddf-ab1c-4fa7-9c51-ba187d647618"
 Inti = "fb74042b-437e-4c5b-88cf-d4e2beb394d5"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
diff --git a/docs/make.jl b/docs/make.jl
index f1710cad..45ad7f80 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,38 +1,56 @@
 using Inti
 using Documenter
-
-# load package needed for extensions
+using Literate
+# packages needed for extensions
 using Gmsh
 using WriteVTK
 using CairoMakie
 using HMatrices
 
-DocMeta.setdocmeta!(Inti, :DocTestSetup, :(using Inti); recursive=true)
+# Generate examples using Literate
+const examples_dir = joinpath(Inti.PROJECT_ROOT, "docs", "src", "examples")
+const generated_dir = joinpath(Inti.PROJECT_ROOT, "docs", "src", "examples", "generated")
+example = "mock_example.jl"
+for example in ["mock_example.jl"]
+    src = joinpath(examples_dir, example)
+    Literate.markdown(src, generated_dir; mdstrings = true)
+    Literate.notebook(src, generated_dir; mdstrings = true)
+end
+
+## setup documentation config
+DocMeta.setdocmeta!(Inti, :DocTestSetup, :(using Inti); recursive = true)
+
+modules = [Inti]
+for extension in [:IntiGmshExt, :IntiMakieExt, :IntiVTKExt, :IntiHMatricesExt]
+    ext = Base.get_extension(Inti, extension)
+    isnothing(ext) && "error loading $ext"
+    push!(modules, ext)
+end
+
+# some settings are only active on a CI build
+on_CI = get(ENV, "CI", "false") == "true"
 
 makedocs(;
-    modules=[
-        Inti,
-    ],
-    authors="Luiz M. Faria",
-    repo="",
-    sitename="Inti.jl",
-    format=Documenter.HTML(;
-        prettyurls=get(ENV, "CI", "false") == "true",
-        canonical="https://IntegralEquations.github.io/Inti.jl",
-        edit_link="main",
+    modules = modules,
+    repo = "",
+    sitename = "Inti.jl",
+    format = Documenter.HTML(;
+        prettyurls = on_CI,
+        canonical = "https://IntegralEquations.github.io/Inti.jl",
+        edit_link = "main",
     ),
-    pages=[
+    pages = [
         "Home" => "index.md",
         "Meshing" => "geo_and_meshes.md",
-        "Examples" =>[
-            "examples/helmholtz_soundsoft_scattering_circle.md"
-        ],
+        "Examples" => ["examples/generated/mock_example.md"],
         "References" => "references.md",
     ],
+    warnonly = on_CI ? false : Documenter.except(:linkcheck_remotes),
     pagesonly = true,
 )
 
 deploydocs(;
-    repo="github.com/IntegralEquations/Inti.jl",
-    devbranch="main",
+    repo = "github.com/IntegralEquations/Inti.jl",
+    devbranch = "main",
+    push_preview = false,
 )
diff --git a/docs/src/examples/helmholtz_soundsoft_scattering_circle.md b/docs/src/examples/helmholtz_soundsoft_scattering_circle.md
deleted file mode 100644
index eebf545e..00000000
--- a/docs/src/examples/helmholtz_soundsoft_scattering_circle.md
+++ /dev/null
@@ -1,206 +0,0 @@
-# Helmholtz sound-soft scattering
-
-```@meta
-CurrentModule = Inti
-```
-
-!!! note "Important points covered in this tutorial"
-    - Creating a geometry and a mesh using `Gmsh`'s API
-    - Assembling integral operators and integral potential
-    - Solving a boundary integral equation
-    - Visualizing the solution
-
-In this tutorial we will show how to solve an acoustic scattering problem in the
-context of Helmholtz equation. We will begin with a simple example involving a
-*smooth* sound-soft obstacle, and gradually make the problem more complex
-considering a transmission problem. Along the way we will introduce the necessary
-techniques used to handle some difficulties encountered. Assuming you
-have already installed `Inti`, let us begin by importing the necessary
-packages for this tutorial:
-
-```@example helmholtz_scattering_2d
-# Load required packages
-using Inti
-using LinearAlgebra
-using StaticArrays
-using CairoMakie
-using Gmsh
-```
-
-We are now ready to begin solving some PDEs!
-
-## Sound-soft scattering
-
-The first example is that of a sound-soft acoustic scattering problem.
-Mathematically, this means we will consider an exterior Helmholtz equation
-(time-harmonic acoustics) with a Dirichlet boundary condition. More precisely,
-let ``\Omega \subset \mathbb{R}^2`` be a bounded domain, and denote by ``\Gamma
-= \partial \Omega`` its boundary. Then we wish to solve
-
-```math
-    \Delta u + k^2 u = 0 \quad \text{on} \quad \mathbb{R}^2 \setminus \bar{\Omega},
-```
-
-subject to Dirichlet boundary conditions on ``\Gamma``:
-
-```math
-    u(\boldsymbol{x}) = g(\boldsymbol{x}) \quad \text{for} \quad \boldsymbol{x} \in \Gamma.
-```
-
-Here ``g`` is the boundary datum, and ``k`` is a constant (the wavenumber).
-
-For concreteness, we will take ``\Gamma`` to be a disk, and focus on
-the *plane-wave scattering* problem. This means we will seek a solution ``u`` of
-the form ``u = u_s + u_i``, where ``u_i`` is a known incident field, and ``u_s``
-is the scattered field we wish to compute. Using the theory of boundary integral
-equations, we can express ``u_s`` as
-
-```math
-    u_s(\boldsymbol{r}) = \mathcal{D}[\sigma](\boldsymbol{r}) - i k \mathcal{S}[\sigma](\boldsymbol{r}),
-```
-
-where ``\mathcal{S}`` is the so-called single layer potential, ``\mathcal{D}``
-is the double-layer potential, and ``\sigma : \Gamma \to \mathbb{C}`` is a
-surface density. This is an indirect formulation (because ``\sigma`` is an
-*auxiliary* density, not necessarily physical) commonly referred to as a
-*combined field formulation*. Taking the limit ``\mathbb{R}^2 \setminus \bar
-\Omega \ni x \to \Gamma``, it can be shown that the following equation holds on
-``\Gamma``:
-
-```math
-    \left( \frac{\mathrm{I}}{2} + \mathrm{D} - i k \mathrm{S} \right)[\sigma] = g,
-```
-
-where $\mathrm{I}$ is the identity operator, and $\mathrm{S}$ and $\mathrm{D}$ are the single- and double-layer operators. This is the **combined field integral equation** that we will solve next.
-
-As for the boundary data ``g``, using the sound-soft condition (i.e. ``u=0`` on
-the scatterer), it follows that ``u_s = -u_i`` on ``\Gamma``. We are now in a
-position to solve the problem! Let us begin by creating the **geometry** and its
-corresponding **mesh**:
-
-```@example helmholtz_scattering_2d
-h = 0.1
-r = 1
-center = SVector(0,0)
-gmsh.initialize()
-gmsh.option.setNumber("Mesh.MeshSizeMax", h)
-# set gmsh verbosity to 2
-gmsh.option.setNumber("General.Verbosity", 2)
-gmsh.model.mesh.setOrder(2)
-Inti.clear_entities!()
-gmsh.model.occ.addDisk(center[1], center[2], 0, r, r)
-gmsh.model.occ.synchronize()
-gmsh.model.mesh.generate(2)
-Ω = Inti.gmsh_import_domain(; dim=2)
-msh = Inti.gmsh_import_mesh(Ω; dim=2)
-gmsh.finalize()
-```
-
-Note that `msh` contains the mesh of both `Ω` and its boundary `Γ`.
-
-!!! tip "Geometry creation"
-    For simple shapes, it is convenient to directly use the `gmsh` API, as
-    illustrated above. As things get more complex, however, it is
-    preferable to use the `gmsh` UI to create your shapes, save a `.geo` or
-    `.msh` file, and then import the file using [`Inti.gmsh_read_msh`](@ref).
-
-Because we will be employing a Nyström method, we must create a quadrature for
-$\Gamma$. This is done by calling the [`Quadrature`](@ref)
-constructor on a mesh object:
-
-```@example helmholtz_scattering_2d
-Γ = Inti.external_boundary(Ω)
-Q = Inti.Quadrature(view(msh,Γ); qorder=5)
-```
-
-With the `Quadrature` constructed, we now can define the integral operators
-``\mathrm{S}`` and ``\mathrm{D}`` associated with the *Helmholtz* equation:
-
-```@example helmholtz_scattering_2d
-k = 2π
-pde = Inti.Helmholtz(;dim=2,k)
-G   = Inti.SingleLayerKernel(pde)
-dG  = Inti.DoubleLayerKernel(pde)
-Sop = Inti.IntegralOperator(G,Q,Q)
-Dop = Inti.IntegralOperator(dG,Q,Q)
-```
-
-Both `Sop` and `Dop` are of type [`IntegralOperator`](@ref), which is a subtype
-of `AbstractMatrix`. They represent a discrete approximation to linear operators
-mapping densities defined on source surface into densities defined on a target
-surface. There are two well-known difficulties related to the discretization of
-these `IntegralOperator`s:
-
-- The kernel of the integral operator is not smooth, and thus specialized
-quadrature rules are required to accurately approximate the matrix entries for
-which the target and source point lie *close* (relative to some scale) to each
-other.
-- The underlying matrix is dense, and thus the storage and computational cost of
-the operator is prohibitive for large problems unless acceleration techniques
-such as *Fast Multipole Methods* or *Hierarchical Matrices* are employed.  
-
-In this example, we will use a sparse correction method based on the
-*density interpolation technique*, and ignore the second issue. More precisely,
-because the problem is two-dimensional and simple, we will just assemble a dense
-`Matrix` to represent the integral operator.
-
-We first build the *dense* part of the operators:
-
-```@example helmholtz_scattering_2d
-S₀ = Matrix(Sop)
-D₀ = Matrix(Dop)
-```
-
-Next we build the sparse corrections
-
-```@example helmholtz_scattering_2d
-δS, δD = Inti.bdim_correction(pde,Q,Q,S₀,D₀)
-```
-
-We can now add the corrections to the dense part to obtain the final operators:
-
-```@example helmholtz_scattering_2d
-S = S₀ + δS
-D = D₀ + δD
-```
-
-Finally, we are ready to solve the scattering problem:
-
-```@example helmholtz_scattering_2d
-# the linear operator
-L = 0.5*I + D - im*k*S
-# incident wave
-θ  = 0*π/4
-d⃗  = (cos(θ),sin(θ))
-uᵢ  = x -> exp(im*k*dot(d⃗,x))
-rhs = [-uᵢ(q.coords) for q in Q]
-σ = L \ rhs
-```
-
-We can now reconstruct the solution using the `IntegralPotential` representation:
-
-```@example helmholtz_scattering_2d
-𝒮 = Inti.IntegralPotential(G,Q)
-𝒟 = Inti.IntegralPotential(dG,Q)
-u = x -> 𝒟[σ](x) - im*k*𝒮[σ](x) + uᵢ(x)
-```
-
-Note that, albeit not very efficient, the function `u` can be evaluated at any
-point. To see what the solution looks like, let's plot it on a grid
-
-```@example helmholtz_scattering_2d
-fig = Figure()
-ax  = Axis(fig[1,1];xgridvisible=false,ygridvisible=false)
-h = 2π/k/10 # 20 pts per wavelength
-x = y = -5:h:5
-# U = [Inti.isinside((x,y),Q) ? NaN + NaN*im : uᵢ((x,y)) for x in x, y in y]
-U = [norm(SVector(x,y) - center) < r ? NaN + NaN*im : u((x,y)) for x in x, y in y]
-hm = heatmap!(ax,x,y,real(U),interpolate=true,colorrange=(-1.0,1.0))
-cm = Colorbar(fig[1, 2],hm; label="Re(u)")
-# plot the geometry
-lines!(ax,view(msh,Γ),linewidth=5.0,color=:black)
-# some final tweak before displaying the figure
-colsize!(fig.layout, 1, Aspect(1, 1.0))
-resize_to_layout!(fig)
-fig
-```
diff --git a/docs/src/examples/mock_example.jl b/docs/src/examples/mock_example.jl
new file mode 100644
index 00000000..5d631143
--- /dev/null
+++ b/docs/src/examples/mock_example.jl
@@ -0,0 +1,22 @@
+import Pkg               #src
+docsdir = joinpath(@__DIR__,"../..") #src
+Pkg.activate(docsdir) #src
+
+md"""
+# Mock example
+"""
+
+#md # [![ipynb](https://img.shields.io/badge/download-ipynb-blue)](mock_example.ipynb)
+#md # [![nbviewer](https://img.shields.io/badge/show-nbviewer-blue.svg)](@__NBVIEWER_ROOT_URL__/examples/generated/mock_example.ipynb)
+
+md"""
+
+Testing how to use `Literate.jl` with `Documenter.jl` to generate documentation.
+
+!!! note "Mock admonition"
+    - This should work using `Documenter Markdown` syntax
+    - If it does not, check back again
+
+"""
+
+using Inti