diff --git a/demos/gray_scott.py b/demos/gray_scott.py
index 2edc961c..36dd3e12 100644
--- a/demos/gray_scott.py
+++ b/demos/gray_scott.py
@@ -159,18 +159,12 @@ def qoi():
 
 # Finally, plot the outputs to be viewed in Paraview. ::
 
-ic = mesh_seq.get_initial_condition()
-for field, sols in solutions.items():
-    fwd_outfile = VTKFile(f"gray_scott/{field}_forward.pvd")
-    adj_outfile = VTKFile(f"gray_scott/{field}_adjoint.pvd")
-    fwd_outfile.write(*ic[field].subfunctions)
-    for i in range(num_subintervals):
-        for sol in sols["forward"][i]:
-            fwd_outfile.write(*sol.subfunctions)
-        for sol in sols["adjoint"][i]:
-            adj_outfile.write(*sol.subfunctions)
-    adj_end = Function(ic[field]).assign(0.0)
-    adj_outfile.write(*adj_end.subfunctions)
+solutions.export(
+    "gray_scott/forward.pvd",
+    export_field_types="forward",
+    initial_condition=mesh_seq.get_initial_condition(),
+)
+solutions.export("gray_scott/adjoint.pvd", export_field_types="adjoint")
 
 # In the `next demo <./gray_scott_split.py.html>`__, we consider solving the same
 # problem, but splitting the solution field into multiple components.
diff --git a/demos/gray_scott_split.py b/demos/gray_scott_split.py
index a44151b0..f8f74b63 100644
--- a/demos/gray_scott_split.py
+++ b/demos/gray_scott_split.py
@@ -161,16 +161,11 @@ def qoi():
 )
 solutions = mesh_seq.solve_adjoint()
 
-ic = mesh_seq.get_initial_condition()
-for field, sols in solutions.items():
-    fwd_outfile = VTKFile(f"gray_scott_split/{field}_forward.pvd")
-    adj_outfile = VTKFile(f"gray_scott_split/{field}_adjoint.pvd")
-    fwd_outfile.write(ic[field])
-    for i in range(num_subintervals):
-        for sol in sols["forward"][i]:
-            fwd_outfile.write(sol)
-        for sol in sols["adjoint"][i]:
-            adj_outfile.write(sol)
-    adj_outfile.write(Function(ic[field]).assign(0.0))
+solutions.export(
+    "gray_scott_split/forward.pvd",
+    export_field_types="forward",
+    initial_condition=mesh_seq.get_initial_condition(),
+)
+solutions.export("gray_scott_split/adjoint.pvd", export_field_types="adjoint")
 
 # This tutorial can be dowloaded as a `Python script <gray_scott_split.py>`__.
diff --git a/demos/solid_body_rotation.py b/demos/solid_body_rotation.py
index acf590aa..b95c0cdb 100644
--- a/demos/solid_body_rotation.py
+++ b/demos/solid_body_rotation.py
@@ -237,14 +237,9 @@ def qoi():
 # use Paraview. To save all adjoint solution components in Paraview format, use the
 # following. ::
 
-for field, sols in solutions.items():
-    fwd_outfile = VTKFile(f"solid_body_rotation/{field}_forward.pvd")
-    adj_outfile = VTKFile(f"solid_body_rotation/{field}_adjoint.pvd")
-    for i in range(len(mesh_seq)):
-        for sol in sols["forward"][i]:
-            fwd_outfile.write(sol)
-        for sol in sols["adjoint"][i]:
-            adj_outfile.write(sol)
+solutions.export(
+    "solid_body_rotation/solutions.pvd", export_field_types=["forward", "adjoint"]
+)
 
 # In the `next demo <./gray_scott.py.html>`__, we increase the complexity by considering
 # two concentration fields in an advection-diffusion-reaction problem.
diff --git a/test_adjoint/test_demos.py b/test_adjoint/test_demos.py
index b46f4ca6..da9a58e1 100644
--- a/test_adjoint/test_demos.py
+++ b/test_adjoint/test_demos.py
@@ -25,16 +25,16 @@
     "burgers-hessian.py": {""""maxiter": 35""": """"maxiter": 3"""},
     "gray_scott.py": {
         "end_time = 2000.0": "end_time = 10.0",
-        r"ic = mesh_seq.get_initial_condition\(\)\n.*?adj_outfile.write\(\*adj_end.subfunctions\)": "",
+        r"solutions\.export\((.*?)\)\s*solutions\.export\((.*?)\)": "",
     },
     "gray_scott_split.py": {
         "end_time = 2000.0": "end_time = 10.0",
-        r"ic = mesh_seq.get_initial_condition\(\)\n.*?adj_outfile.write\(Function\(ic\[field\]\).assign\(0.0\)\)": "",
+        r"solutions\.export\((.*?)\)\s*solutions\.export\((.*?)\)": "",
     },
     "point_discharge2d-hessian.py": {""""maxiter": 35""": """"maxiter": 3"""},
     "point_discharge2d-goal_oriented.py": {""""maxiter": 35""": """"maxiter": 3"""},
     "solid_body_rotation.py": {
-        r"for field, sols in solutions.items\(\):\n.*?adj_outfile.write\(sol\)": "",
+        r"solutions\.export\((.*?)\)": "",
     },
 }