diff --git a/cpp/demo/CMakeLists.txt b/cpp/demo/CMakeLists.txt
index dd4cf14d928..6c3b28524f6 100644
--- a/cpp/demo/CMakeLists.txt
+++ b/cpp/demo/CMakeLists.txt
@@ -20,3 +20,5 @@ add_demo_subdirectory(poisson)
 add_demo_subdirectory(poisson_matrix_free)
 add_demo_subdirectory(hyperelasticity)
 add_demo_subdirectory(interpolation-io)
+add_demo_subdirectory(interpolation_different_meshes)
+
diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
index 377f9d78f5e..5b2bf8be35a 100644
--- a/cpp/demo/interpolation_different_meshes/main.cpp
+++ b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -5,13 +5,13 @@
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
 #include <basix/e-lagrange.h>
-#include <dolfinx/common/timing.h>
-#include <dolfinx/fem/utils.h>
+#include <dolfinx/fem/dolfin_fem.h>
 #include <dolfinx/io/ADIOS2Writers.h>
 #include <dolfinx/mesh/generation.h>
 #include <memory>
 
 using namespace dolfinx;
+namespace stdex = std::experimental;
 
 using T = double;
 
@@ -24,38 +24,46 @@ int main(int argc, char* argv[])
 
     // Create a tetrahedral mesh
     auto mesh_tet = std::make_shared<mesh::Mesh>(
-        mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {10, 10, 10},
-                         mesh::CellType::tetrahedron, mesh::GhostMode::none));
+        mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {20, 20, 20},
+                         mesh::CellType::tetrahedron));
     // Create a hexahedral mesh
     auto mesh_hex = std::make_shared<mesh::Mesh>(
-        mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {9, 8, 7},
-                         mesh::CellType::hexahedron, mesh::GhostMode::none));
+        mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {15, 15, 15},
+                         mesh::CellType::hexahedron));
 
-    basix::FiniteElement eL = basix::element::create_lagrange(
+    basix::FiniteElement element_tet = basix::element::create_lagrange(
         mesh::cell_type_to_basix_type(mesh_tet->topology().cell_type()), 1,
         basix::element::lagrange_variant::equispaced, false);
     auto V_tet = std::make_shared<fem::FunctionSpace>(
-        fem::create_functionspace(mesh_tet, eL, 3));
+        fem::create_functionspace(mesh_tet, element_tet, 3));
 
-    basix::FiniteElement eR = basix::element::create_lagrange(
+    basix::FiniteElement element_hex = basix::element::create_lagrange(
         mesh::cell_type_to_basix_type(mesh_hex->topology().cell_type()), 2,
         basix::element::lagrange_variant::equispaced, false);
     auto V_hex = std::make_shared<fem::FunctionSpace>(
-        fem::create_functionspace(mesh_hex, eR, 3));
+        fem::create_functionspace(mesh_hex, element_hex, 3));
 
     auto u_tet = std::make_shared<fem::Function<T>>(V_tet);
     auto u_hex = std::make_shared<fem::Function<T>>(V_hex);
 
-    auto fun = [](auto& x)
+    auto fun = [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
     {
-      auto r = xt::zeros_like(x);
-      xt::row(r, 0) = xt::cos(10 * xt::row(x, 0)) * xt::sin(10 * xt::row(x, 2));
-      xt::row(r, 1) = xt::sin(10 * xt::row(x, 0)) * xt::sin(10 * xt::row(x, 2));
-      xt::row(r, 2) = xt::cos(10 * xt::row(x, 0)) * xt::cos(10 * xt::row(x, 2));
-      return r;
+      std::vector<T> fdata(3 * x.extent(1), 0.0);
+      using dextent = stdex::dextents<std::size_t, 2>;
+      stdex::mdspan<double, dextent> f(fdata.data(), 3, x.extent(1));
+      for (std::size_t i = 0; i < x.extent(1); ++i)
+      {
+        f(0, i) = std::cos(10 * x(0, i)) * std::sin(10 * x(2, i));
+        f(1, i) = std::sin(10 * x(0, i)) * std::sin(10 * x(2, i));
+        f(2, i) = std::cos(10 * x(0, i)) * std::cos(10 * x(2, i));
+      }
+      return {std::move(fdata), {3, x.extent(1)}};
     };
 
+    // Interpolate an expression into u_tet
     u_tet->interpolate(fun);
+
+    // Interpolate from u_tet to u_hex
     u_hex->interpolate(*u_tet);
 
 #ifdef HAS_ADIOS2