diff --git a/doc/userdoc.org b/doc/userdoc.org
index da3ac7c0192a1c3fb214a7811a7525bc95e5c888..9295c530bbf2773a6fcc2fe67e1ec9d07c52fd15 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -3218,11 +3218,6 @@ Here is the list of the functions
 - ~tan: Vh -> Vh~
 - ~tanh: Vh -> Vh~
 
-Additionally the following function is defined from
-$\vec{\mathbb{P}}_0(\mathbb{R})$ to $\mathbb{P}_0(\mathbb{R})$
-- ~sum_to_Vh: Vh -> Vh~
-it computes the sum of all the components of the vector function.
-
 ******  ~Vh*Vh -> Vh~
 
 These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
@@ -3411,6 +3406,19 @@ and the return value is an $\mathbb{R}^{3\times3}$ matrix.
   return value is $$\sum_{j\in\mathcal{J}} A_j.$$
 
 
+***** Utilities
+
+****** ~sum_to_Vh: Vh -> Vh~
+
+This function it computes the sum of all the components of
+$\vec{\mathbb{P}}_0(\mathbb{R})$ function and stores the result into a
+$\mathbb{P}_0(\mathbb{R})$.
+
+****** ~vectorize: (Vh) -> Vh~
+
+This function creates a $\vec{\mathbb{P}}_0(\mathbb{R})$ function using
+a tuple of $\mathbb{P}_0(\mathbb{R})$.
+
 ***** Interpolation and integration functions
 
 These functions are ways to define discrete functions from analytic
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index 99d9bfcec67af3819bf38236b9f0348149f80baa..fcf3b31a73da2c47bb59e27bc169195255effb30 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -331,6 +331,15 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
 
                                       ));
 
+  scheme_module
+    ._addBuiltinFunction("vectorize",
+                         std::function(
+
+                           [](const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list)
+                             -> std::shared_ptr<const IDiscreteFunction> { return vectorize(discrete_function_list); }
+
+                           ));
+
   scheme_module._addBuiltinFunction("integral_of_R", std::function(
 
                                                        [](std::shared_ptr<const IDiscreteFunction> a) -> double {
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index 25ec92c298c21052e61003d13fe1edaa13106ce1..f83bd418592c29b5cded8b77bda6c717ff3af804 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -810,6 +810,91 @@ sum_of_Vh_components(const std::shared_ptr<const IDiscreteFunction>& f)
   }
 }
 
+std::shared_ptr<const IDiscreteFunction>
+vectorize(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list)
+{
+  Assert(discrete_function_list.size() > 0);
+  std::shared_ptr p_i_mesh = [&] {
+    auto i = discrete_function_list.begin();
+
+    std::shared_ptr<const IMesh> i_mesh = (*i)->mesh();
+    ++i;
+    for (; i != discrete_function_list.end(); ++i) {
+      if ((*i)->mesh() != i_mesh) {
+        throw NormalError("discrete functions are not defined on the same mesh");
+      }
+    }
+
+    return i_mesh;
+  }();
+
+  for (auto&& i_discrete_function : discrete_function_list) {
+    if ((i_discrete_function->descriptor().type() != DiscreteFunctionType::P0) or
+        (i_discrete_function->dataType() != ASTNodeDataType::double_t)) {
+      throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(i_discrete_function));
+    }
+  }
+
+  switch (p_i_mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension       = 1;
+    using DiscreteFunctionVectorType = DiscreteFunctionP0Vector<Dimension, double>;
+    using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, double>;
+    std::shared_ptr<const Mesh<Connectivity<Dimension>>> p_mesh =
+      std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(p_i_mesh);
+
+    DiscreteFunctionVectorType vector_function(p_mesh, discrete_function_list.size());
+    for (size_t i = 0; i < discrete_function_list.size(); ++i) {
+      const DiscreteFunctionType& f = dynamic_cast<const DiscreteFunctionType&>(*discrete_function_list[i]);
+      for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        vector_function[cell_id][i] = f[cell_id];
+      }
+    }
+
+    return std::make_shared<DiscreteFunctionVectorType>(vector_function);
+  }
+  case 2: {
+    constexpr size_t Dimension       = 2;
+    using DiscreteFunctionVectorType = DiscreteFunctionP0Vector<Dimension, double>;
+    using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, double>;
+    std::shared_ptr<const Mesh<Connectivity<Dimension>>> p_mesh =
+      std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(p_i_mesh);
+
+    DiscreteFunctionVectorType vector_function(p_mesh, discrete_function_list.size());
+    for (size_t i = 0; i < discrete_function_list.size(); ++i) {
+      const DiscreteFunctionType& f = dynamic_cast<const DiscreteFunctionType&>(*discrete_function_list[i]);
+      for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        vector_function[cell_id][i] = f[cell_id];
+      }
+    }
+
+    return std::make_shared<DiscreteFunctionVectorType>(vector_function);
+  }
+  case 3: {
+    constexpr size_t Dimension       = 3;
+    using DiscreteFunctionVectorType = DiscreteFunctionP0Vector<Dimension, double>;
+    using DiscreteFunctionType       = DiscreteFunctionP0<Dimension, double>;
+    std::shared_ptr<const Mesh<Connectivity<Dimension>>> p_mesh =
+      std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(p_i_mesh);
+
+    DiscreteFunctionVectorType vector_function(p_mesh, discrete_function_list.size());
+    for (size_t i = 0; i < discrete_function_list.size(); ++i) {
+      const DiscreteFunctionType& f = dynamic_cast<const DiscreteFunctionType&>(*discrete_function_list[i]);
+      for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        vector_function[cell_id][i] = f[cell_id];
+      }
+    }
+
+    return std::make_shared<DiscreteFunctionVectorType>(vector_function);
+  }
+    // LCOV_EXCL_START
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+    // LCOV_EXCL_STOP
+  }
+}
+
 template double sum_of<double>(const std::shared_ptr<const IDiscreteFunction>&);
 
 template TinyVector<1> sum_of<TinyVector<1>>(const std::shared_ptr<const IDiscreteFunction>&);
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
index 360175e3ca416f6e0058b77076748823759bb346..2369429a6267a8fdbc803d75b4213703dd970a48 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -67,6 +67,9 @@ std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&,
 
 std::shared_ptr<const IDiscreteFunction> sum_of_Vh_components(const std::shared_ptr<const IDiscreteFunction>&);
 
+std::shared_ptr<const IDiscreteFunction> vectorize(
+  const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list);
+
 double min(const std::shared_ptr<const IDiscreteFunction>&);
 
 std::shared_ptr<const IDiscreteFunction> min(const std::shared_ptr<const IDiscreteFunction>&,
diff --git a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
index 5370595face3f88f7dfbffdfe186c53d3b5aa612..5386d9ae69f9b7b4fbd093d82e1daaef2cbc81c7 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -571,6 +571,35 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(sum_of_Vh_components(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
         }
 
+        SECTION("vectorize (Vh) -> Vh")
+        {
+          {
+            std::shared_ptr p_sum_components =
+              vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_positive_u, p_bounded_u});
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0Vector<Dimension, double> sum_components =
+              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_sum_components);
+            REQUIRE(sum_components.size() == 3);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              is_same &= ((*p_u)[cell_id] == sum_components[cell_id][0]);
+              is_same &= ((*p_positive_u)[cell_id] == sum_components[cell_id][1]);
+              is_same &= ((*p_bounded_u)[cell_id] == sum_components[cell_id][2]);
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_other_mesh_u}),
+                              "error: discrete functions are not defined on the same mesh");
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_R1_u}),
+                              "error: invalid operand type Vh(P0:R^1)");
+        }
+
         SECTION("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
@@ -1114,6 +1143,33 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(sum_of_Vh_components(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
         }
 
+        SECTION("vectorize (Vh) -> Vh")
+        {
+          {
+            std::shared_ptr p_sum_components =
+              vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_positive_u, p_bounded_u});
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0Vector<Dimension, double> sum_components =
+              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_sum_components);
+            REQUIRE(sum_components.size() == 3);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              is_same &= ((*p_u)[cell_id] == sum_components[cell_id][0]);
+              is_same &= ((*p_positive_u)[cell_id] == sum_components[cell_id][1]);
+              is_same &= ((*p_bounded_u)[cell_id] == sum_components[cell_id][2]);
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_other_mesh_u}),
+                              "error: discrete functions are not defined on the same mesh");
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_R1_u}),
+                              "error: invalid operand type Vh(P0:R^1)");
+        }
+
         SECTION("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);
@@ -1657,6 +1713,35 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           REQUIRE_THROWS_WITH(sum_of_Vh_components(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
         }
 
+        SECTION("vectorize (Vh) -> Vh")
+        {
+          {
+            std::shared_ptr p_sum_components =
+              vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_positive_u, p_bounded_u});
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            REQUIRE(p_sum_components.use_count() == 1);
+
+            const DiscreteFunctionP0Vector<Dimension, double> sum_components =
+              dynamic_cast<const DiscreteFunctionP0Vector<Dimension, double>&>(*p_sum_components);
+            REQUIRE(sum_components.size() == 3);
+
+            DiscreteFunctionP0<Dimension, double> direct_sum(mesh);
+            bool is_same = true;
+            for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+              is_same &= ((*p_u)[cell_id] == sum_components[cell_id][0]);
+              is_same &= ((*p_positive_u)[cell_id] == sum_components[cell_id][1]);
+              is_same &= ((*p_bounded_u)[cell_id] == sum_components[cell_id][2]);
+            }
+            REQUIRE(is_same);
+          }
+
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_u, p_other_mesh_u}),
+                              "error: discrete functions are not defined on the same mesh");
+          REQUIRE_THROWS_WITH(vectorize(std::vector<std::shared_ptr<const IDiscreteFunction>>{p_R1_u}),
+                              "error: invalid operand type Vh(P0:R^1)");
+        }
+
         SECTION("dot Vh*Rd -> Vh")
         {
           CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), dot);