diff --git a/docs/source/interpolation.rst b/docs/source/interpolation.rst
index 20a4467714..56d972aac9 100644
--- a/docs/source/interpolation.rst
+++ b/docs/source/interpolation.rst
@@ -151,7 +151,14 @@ function space:
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 11-23, 31-40
+   :start-after: [test_line_integral 1]
+   :end-before: [test_line_integral 2]
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_line_integral 3]
+   :end-before: [test_line_integral 4]
 
 For more on forms, see :ref:`this section of the manual <more_complicated_forms>`.
 
@@ -170,7 +177,14 @@ interpolation will raise a :py:class:`~.DofNotDefinedError`.
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 46-58, 64-65
+   :start-after: [test_cross_mesh 1]
+   :end-before: [test_cross_mesh 2]
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_cross_mesh 3]
+   :end-before: [test_cross_mesh 4]
 
 This can be overriden with the optional ``allow_missing_dofs`` keyword
 argument:
@@ -178,7 +192,14 @@ argument:
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 72-73, 80
+   :start-after: [test_cross_mesh 5]
+   :end-before: [test_cross_mesh 6]
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_cross_mesh 7]
+   :end-before: [test_cross_mesh 8]
 
 In this case, the missing degrees of freedom (DoFs, the global basis function
 coefficients which could not be set) are, by default, set to zero:
@@ -186,7 +207,8 @@ coefficients which could not be set) are, by default, set to zero:
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 85
+   :start-after: [test_cross_mesh 9]
+   :end-before: [test_cross_mesh 10]
 
 If we specify an output :py:class:`~.Function` then the missing DoFs are
 unmodified.
@@ -197,7 +219,8 @@ we set them to be ``nan`` ('not a number') for easy identification:
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 90-93
+   :start-after: [test_cross_mesh 11]
+   :end-before: [test_cross_mesh 12]
 
 If we specify an output :py:class:`~.Function`, this overwrites the missing
 DoFs.
@@ -208,7 +231,8 @@ argument is set at construction:
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 100
+   :start-after: [test_cross_mesh 13]
+   :end-before: [test_cross_mesh 14]
 
 The ``default_missing_val`` keyword argument is then set whenever we call
 :py:meth:`~.Interpolator.interpolate`:
@@ -216,7 +240,8 @@ The ``default_missing_val`` keyword argument is then set whenever we call
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 103
+   :start-after: [test_cross_mesh 15]
+   :end-before: [test_cross_mesh 16]
 
 If we supply an output :py:class:`~.Function` and don't set
 ``default_missing_val`` then any missing DoFs are left as they were prior to
@@ -225,7 +250,20 @@ interpolation:
 .. literalinclude:: ../../tests/regression/test_interpolation_manual.py
    :language: python3
    :dedent:
-   :lines: 107-110, 113-115, 124-126
+   :start-after: [test_cross_mesh 17]
+   :end-before: [test_cross_mesh 18]
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_cross_mesh 19]
+   :end-before: [test_cross_mesh 20]
+
+.. literalinclude:: ../../tests/regression/test_interpolation_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_cross_mesh 21]
+   :end-before: [test_cross_mesh 22]
 
 
 Interpolation from external data
diff --git a/docs/source/point-evaluation.rst b/docs/source/point-evaluation.rst
index 9357dc140f..b14b459c69 100644
--- a/docs/source/point-evaluation.rst
+++ b/docs/source/point-evaluation.rst
@@ -131,7 +131,8 @@ evaluation of a function :math:`f` defined in a function space
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 9-28
+   :start-after: [test_vertex_only_mesh_manual_example 1]
+   :end-before: [test_vertex_only_mesh_manual_example 3]
 
 will print ``[0.02, 0.08, 0.18]`` when running in serial, the values of
 :math:`x^2 + y^2` at the points :math:`(0.1, 0.1)`, :math:`(0.2, 0.2)` and
@@ -225,7 +226,20 @@ warning or switched off entirely:
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 50-52,55-57,59-63
+   :start-after: [test_vom_manual_points_outside_domain 1]
+   :end-before: [test_vom_manual_points_outside_domain 2]
+
+.. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_vom_manual_points_outside_domain 3]
+   :end-before: [test_vom_manual_points_outside_domain 4]
+
+.. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_vom_manual_points_outside_domain 5]
+   :end-before: [test_vom_manual_points_outside_domain 6]
 
 
 Expressions with point evaluations
@@ -254,7 +268,8 @@ Firedrake using :func:`~.VertexOnlyMesh` and :func:`~.interpolate` as
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 72-86
+   :start-after: [test_vom_manual_keyword_arguments 1]
+   :end-before: [test_vom_manual_keyword_arguments 2]
 
 .. _external-point-data:
 
@@ -276,7 +291,8 @@ vertex-only mesh. For example:
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 158-172
+   :start-after: [test_input_ordering_input 1]
+   :end-before: [test_input_ordering_input 2]
 
 This is entirely parallel safe.
 
@@ -288,7 +304,8 @@ had
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 26-28
+   :start-after: [test_vertex_only_mesh_manual_example 2]
+   :end-before: [test_vertex_only_mesh_manual_example 3]
 
 In parallel, this will print the values of ``f`` at the given ``points`` list
 **after the points have been distributed over the parent mesh**. If we want the
@@ -298,7 +315,8 @@ distributed** we can use :py:attr:`~.VertexOnlyMeshTopology.input_ordering` as f
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 30-37
+   :start-after: [test_vertex_only_mesh_manual_example 4]
+   :end-before: [test_vertex_only_mesh_manual_example 5]
 
 .. note::
 
@@ -321,7 +339,8 @@ a good idea to set the values to ``nan`` before the interpolation:
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 39-43
+   :start-after: [test_vertex_only_mesh_manual_example 6]
+   :end-before: [test_vertex_only_mesh_manual_example 7]
 
 
 More ways to interact with external data
@@ -363,7 +382,8 @@ the mesh cells and is a property of the mesh itself
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 92-112
+   :start-after: [test_mesh_tolerance 1]
+   :end-before: [test_mesh_tolerance 2]
 
 Keyword arguments
 ~~~~~~~~~~~~~~~~~
@@ -380,7 +400,8 @@ vertex-only mesh. This will modify the tolerance property of the parent mesh.
 .. literalinclude:: ../../tests/vertexonly/test_vertex_only_manual.py
    :language: python3
    :dedent:
-   :lines: 122-139
+   :start-after: [test_mesh_tolerance_change 1]
+   :end-before: [test_mesh_tolerance_change 2]
 
 Note that since our tolerance is relative, the number of cells in a mesh
 dictates the point loss behaviour close to cell edges. So the mesh
diff --git a/tests/regression/test_interpolation_manual.py b/tests/regression/test_interpolation_manual.py
index 348e6e4591..46164b6f3b 100644
--- a/tests/regression/test_interpolation_manual.py
+++ b/tests/regression/test_interpolation_manual.py
@@ -3,12 +3,9 @@
 import pytest
 import numpy as np
 
-# Ensure that the code shown in the manual runs without error. If you change
-# the code here make sure you update the .. literalinclude:: bits in the manual
-# too!
-
 
 def test_line_integral():
+    # [test_line_integral 1]
     # Start with a simple field exactly represented in the function space over
     # the unit square domain.
     m = UnitSquareMesh(2, 2)
@@ -22,6 +19,7 @@ def test_line_integral():
     vertex_coords = np.asarray([[0.0, 0.0], [1.0, 1.0]])
     plex = mesh.plex_from_cell_list(1, cells, vertex_coords, comm=m.comm)
     line = mesh.Mesh(plex, dim=2)
+    # [test_line_integral 2]
     x, y = SpatialCoordinate(line)
     V_line = FunctionSpace(line, "CG", 2)
     f_line = Function(V_line).interpolate(x * y)
@@ -29,6 +27,7 @@ def test_line_integral():
     f_line.zero()
     assert np.isclose(assemble(f_line * dx), 0)  # sanity again
 
+    # [test_line_integral 3]
     # We want to calculate the line integral of f along it. To do this we
     # create a function space on the line mesh...
     V_line = FunctionSpace(line, "CG", 2)
@@ -39,11 +38,13 @@ def test_line_integral():
     # The integral of f along the line is then a simple form expression which
     # we assemble:
     assemble(f_line * dx)  # this outputs sqrt(2) / 3
+    # [test_line_integral 4]
     assert np.isclose(assemble(f_line * dx), np.sqrt(2) / 3)
 
 
 def test_cross_mesh():
     def correct_indent():
+        # [test_cross_mesh 1]
         # These meshes only share some of their domain
         src_mesh = UnitSquareMesh(2, 2)
         dest_mesh = UnitSquareMesh(3, 3, quadrilateral=True)
@@ -56,26 +57,32 @@ def correct_indent():
 
         # ... and want to interpolate into a function space on our target mesh ...
         V_dest = FunctionSpace(dest_mesh, "Q", 2)
-
+        # [test_cross_mesh 2]
         return src_mesh, dest_mesh, f_src, V_dest
 
     src_mesh, dest_mesh, f_src, V_dest = correct_indent()
 
     with pytest.raises(DofNotDefinedError):
+        # [test_cross_mesh 3]
         # ... but get a DofNotDefinedError if we try
         f_dest = assemble(interpolate(f_src, V_dest))  # raises DofNotDefinedError
+        # [test_cross_mesh 4]
 
     with pytest.raises(DofNotDefinedError):
         # as will the interpolate method of a Function
         f_dest = Function(V_dest).interpolate(f_src)
 
+    # [test_cross_mesh 5]
     # Setting the allow_missing_dofs keyword allows the interpolation to proceed.
     f_dest = assemble(interpolate(f_src, V_dest, allow_missing_dofs=True))
+    # [test_cross_mesh 6]
 
     assert np.isclose(f_dest.at(0.5, 0.5), 0.5)
 
     # or
+    # [test_cross_mesh 7]
     f_dest = Function(V_dest).interpolate(f_src, allow_missing_dofs=True)
+    # [test_cross_mesh 8]
 
     # We get values at the points in the destination mesh as we would expect
     f_dest.at(0.5, 0.5)  # returns 0.5**2 + 0.5**2 = 0.5
@@ -83,36 +90,48 @@ def correct_indent():
     assert np.isclose(f_dest.at(0.5, 0.5), 0.5)
 
     # By default the missing points are set to 0.0
+    # [test_cross_mesh 9]
     f_dest.at(1.5, 1.5)  # returns 0.0
+    # [test_cross_mesh 10]
 
     assert np.isclose(f_dest.at(1.5, 1.5), 0.0)
 
     # We can alternatively specify a value to use for missing points:
+    # [test_cross_mesh 11]
     f_dest = assemble(interpolate(
         f_src, V_dest, allow_missing_dofs=True, default_missing_val=np.nan
     ))
     f_dest.at(1.5, 1.5)  # returns np.nan
+    # [test_cross_mesh 12]
 
     assert np.isclose(f_dest.at(0.5, 0.5), 0.5)
     assert np.isnan(f_dest.at(1.5, 1.5))
 
     # When creating interpolators, the allow_missing_dofs keyword argument is
     # set when creating the interpolator, rather than when calling interpoalate.
+    # [test_cross_mesh 13]
     interpolator = Interpolator(f_src, V_dest, allow_missing_dofs=True)
+    # [test_cross_mesh 14]
 
     # A default missing value can be specified when calling interpolate.
+    # [test_cross_mesh 15]
     f_dest = assemble(interpolator.interpolate(default_missing_val=np.nan))
+    # [test_cross_mesh 16]
 
     # If we supply an output function and don't set default_missing_val
     # then any points outside the domain are left as they were.
+    # [test_cross_mesh 17]
     x_dest, y_dest = SpatialCoordinate(dest_mesh)
     f_dest = Function(V_dest).interpolate(x_dest + y_dest)
     f_dest.at(0.5, 0.5)  # returns x_dest + y_dest = 1.0
+    # [test_cross_mesh 18]
 
     assert np.isclose(f_dest.at(0.5, 0.5), 1.0)
 
+    # [test_cross_mesh 19]
     assemble(interpolator.interpolate(), tensor=f_dest)
     f_dest.at(0.5, 0.5)  # now returns x_src^2 + y_src^2 = 0.5
+    # [test_cross_mesh 20]
 
     assert np.isclose(f_dest.at(0.5, 0.5), 0.5)
 
@@ -122,9 +141,11 @@ def correct_indent():
     f_dest.interpolate(x_dest + y_dest)
     assert np.isclose(f_dest.at(0.5, 0.5), 1.0)  # x_dest + y_dest = 1.0
 
+    # [test_cross_mesh 21]
     # Similarly, using the interpolate method on a function will not overwrite
     # the pre-existing values if default_missing_val is not set
     f_dest.interpolate(f_src, allow_missing_dofs=True)
+    # [test_cross_mesh 22]
 
     assert np.isclose(f_dest.at(0.5, 0.5), 0.5)  # x_src^2 + y_src^2 = 0.5
     assert np.isclose(f_dest.at(1.5, 1.5), 3.0)  # x_dest + y_dest = 3.0
diff --git a/tests/vertexonly/test_vertex_only_manual.py b/tests/vertexonly/test_vertex_only_manual.py
index 80cce57527..7578b5c69b 100644
--- a/tests/vertexonly/test_vertex_only_manual.py
+++ b/tests/vertexonly/test_vertex_only_manual.py
@@ -1,12 +1,10 @@
 from firedrake import *
 from firedrake.__future__ import *
 import pytest
-# Ensure that the code shown in the manual runs without error. If you change
-# the code here make sure you update the .. literalinclude:: bits in the manual
-# too!
 
 
 def test_vertex_only_mesh_manual_example():
+    # [test_vertex_only_mesh_manual_example 1]
     parent_mesh = UnitSquareMesh(10, 10)
 
     V = FunctionSpace(parent_mesh, "CG", 2)
@@ -24,10 +22,13 @@ def test_vertex_only_mesh_manual_example():
     P0DG = FunctionSpace(vom, "DG", 0)
 
     # Interpolation performs point evaluation
+    # [test_vertex_only_mesh_manual_example 2]
     f_at_points = assemble(interpolate(f, P0DG))
 
     print(f_at_points.dat.data_ro)
+    # [test_vertex_only_mesh_manual_example 3]
 
+    # [test_vertex_only_mesh_manual_example 4]
     # Make a P0DG function on the input ordering vertex-only mesh again
     P0DG_input_ordering = FunctionSpace(vom.input_ordering, "DG", 0)
     f_at_input_points = Function(P0DG_input_ordering)
@@ -36,32 +37,41 @@ def test_vertex_only_mesh_manual_example():
     f_at_input_points.interpolate(f_at_points)
 
     print(f_at_input_points.dat.data_ro)  # will print the values at the input points
+    # [test_vertex_only_mesh_manual_example 5]
 
+    # [test_vertex_only_mesh_manual_example 6]
     import numpy as np
     f_at_input_points.dat.data_wo[:] = np.nan
     f_at_input_points.interpolate(f_at_points)
 
     print(f_at_input_points.dat.data_ro)  # any points not found will be NaN
+    # [test_vertex_only_mesh_manual_example 7]
 
 
 def test_vom_manual_points_outside_domain():
     for i in [0]:
         parent_mesh = UnitSquareMesh(100, 100, quadrilateral=True)
 
+        # [test_vom_manual_points_outside_domain 1]
         # point (1.1, 1.0) is outside the mesh
         points = [[0.1, 0.1], [0.2, 0.2], [1.1, 1.0]]
+        # [test_vom_manual_points_outside_domain 2]
 
     vom = True  # avoid flake8 unused variable warning
     with pytest.raises(VertexOnlyMeshMissingPointsError):
+        # [test_vom_manual_points_outside_domain 3]
         # This will raise a VertexOnlyMeshMissingPointsError
         vom = VertexOnlyMesh(parent_mesh, points, missing_points_behaviour='error')
+        # [test_vom_manual_points_outside_domain 4]
 
     def display_correct_indent():
+        # [test_vom_manual_points_outside_domain 5]
         # This will generate a warning and the point will be lost
         vom = VertexOnlyMesh(parent_mesh, points, missing_points_behaviour='warn')
 
         # This will cause the point to be silently lost
         vom = VertexOnlyMesh(parent_mesh, points, missing_points_behaviour=None)
+        # [test_vom_manual_points_outside_domain 6]
 
         assert vom  # Just here to shut up flake8 unused variable warning.
 
@@ -70,6 +80,7 @@ def display_correct_indent():
 
 
 def test_vom_manual_keyword_arguments():
+    # [test_vom_manual_keyword_arguments 1]
     omega = UnitSquareMesh(100, 100, quadrilateral=True)
 
     V = FunctionSpace(omega, "CG", 2)
@@ -86,11 +97,13 @@ def test_vom_manual_keyword_arguments():
     # the points
     f_at_points = assemble(interpolate(f, P0DG))
     expr = assemble(f_at_points * dx)
+    # [test_vom_manual_keyword_arguments 2]
 
     assert expr == 0.0
 
 
 def test_mesh_tolerance():
+    # [test_mesh_tolerance 1]
     parent_mesh = UnitSquareMesh(100, 100, quadrilateral=True)
 
     # point (1.1, 1.0) is outside the mesh
@@ -112,6 +125,7 @@ def test_mesh_tolerance():
     # Similarly .at will not generate an error
     V = FunctionSpace(parent_mesh, 'CG', 2)
     Function(V).at((1.1, 1.0))
+    # [test_mesh_tolerance 2]
 
     assert vom
 
@@ -121,6 +135,7 @@ def test_mesh_tolerance_change():
     points = [[0.1, 0.1], [0.2, 0.2], [1.1, 1.0]]
     V = FunctionSpace(parent_mesh, 'CG', 2)
 
+    # [test_mesh_tolerance_change 1]
     # The point (1.1, 1.0) will still be included in the vertex-only mesh
     vom = VertexOnlyMesh(parent_mesh, points, tolerance=30.0)
 
@@ -139,6 +154,7 @@ def test_mesh_tolerance_change():
     except PointNotInDomainError:
         # But the tolerance property has still been changed - this will print 1.0
         print(parent_mesh.tolerance)
+    # [test_mesh_tolerance_change 2]
 
     assert vom
 
@@ -157,6 +173,7 @@ def test_input_ordering_input():
         point_locations_from_elsewhere = np.array([]).reshape(0, 2)
         point_data_values_from_elsewhere = []
 
+    # [test_input_ordering_input 1]
     # We have a set of points with corresponding data from elsewhere which vary
     # from rank to rank
     vom = VertexOnlyMesh(parent_mesh, point_locations_from_elsewhere, redundant=False)
@@ -172,6 +189,7 @@ def test_input_ordering_input():
 
     # Interpolate puts this data onto the original vertex-only mesh
     point_data = assemble(interpolate(point_data_input_ordering, P0DG))
+    # [test_input_ordering_input 2]
 
     assert vom
     if point_data:  # in case point data is None for some reason - apparently it can be in complex mode without an error occuring?