diff --git a/2-solve-heat-equation.py b/2-solve-heat-equation.py
index 4b98fcf..daad669 100644
--- a/2-solve-heat-equation.py
+++ b/2-solve-heat-equation.py
@@ -40,6 +40,7 @@
     apply_lifting,
     set_bc,
 )
+from pathlib import Path
 
 t = 0.0  # Start time
 T = 10.0  # Final time
@@ -58,14 +59,19 @@
 boundary_facets = mesh.locate_entities_boundary(
     domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool)
 )
-bc = fem.dirichletbc(
-    PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V
-)
+bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)
 
 
 xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
 xdmf.write_mesh(domain)
 
+filename = Path("diffusion.vtkhdf")
+backend_args = {"name": "MyGrid"}
+
+io4dolfinx.write_mesh(
+    filename, domain, backend="vtkhdf", mode=io4dolfinx.FileMode.write, backend_args=backend_args
+)
+
 uh = fem.Function(V)
 uh.name = "uh"
 
@@ -87,7 +93,6 @@
 solver.setType(PETSc.KSP.Type.PREONLY)
 solver.getPC().setType(PETSc.PC.Type.LU)
 
-
 for i in range(num_steps):
     t += dt
     print(f"Time step {i + 1}/{num_steps}, Time: {t:.2f}")
@@ -112,6 +117,15 @@
     # Write solution to file
     xdmf.write_function(uh, t)
 
+    io4dolfinx.write_point_data(
+        filename,
+        uh,
+        time=t,
+        backend="vtkhdf",
+        mode=io4dolfinx.FileMode.append,
+        backend_args=backend_args,
+    )
+
 xdmf.close()
 
 A.destroy()