diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/dataset.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/dataset.hpp
index 3108f6307c..5f5ea8ad38 100644
--- a/power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/dataset.hpp
+++ b/power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/dataset.hpp
@@ -287,6 +287,9 @@ template <dataset_type_tag dataset_type_> class Dataset {
     }
     constexpr bool is_dense(Idx const i) const { return is_dense(buffers_[i]); }
     constexpr bool is_dense(Buffer const& buffer) const { return buffer.indptr.empty(); }
+    constexpr bool is_dense() const {
+        return std::ranges::all_of(buffers_, [this](Buffer const& buffer) { return is_dense(buffer); });
+    }
     constexpr bool is_sparse(std::string_view component, bool with_attribute_buffers = false) const {
         Idx const idx = find_component(component, false);
         if (idx == invalid_index) {
@@ -510,6 +513,35 @@ template <dataset_type_tag dataset_type_> class Dataset {
         return result;
     }
 
+    // get slice dataset from batch
+    Dataset get_slice_scenario(Idx begin, Idx end) const
+        requires(!is_indptr_mutable_v<dataset_type>)
+    {
+        assert(begin <= end);
+        assert(0 <= begin);
+        assert(end <= batch_size());
+        assert(is_batch());
+        assert(is_dense());
+
+        // empty slice
+        if (begin == end) {
+            Dataset result{true, 0, dataset_info_.dataset->name, *meta_data_};
+            result.add_buffer("node", 0, 0, nullptr, nullptr);
+            return result;
+        }
+
+        // start with begin
+        Dataset result = get_individual_scenario(begin);
+        Idx const batch_size = end - begin;
+        result.dataset_info_.is_batch = true;
+        result.dataset_info_.batch_size = batch_size;
+        for (auto&& [buffer, component_info] : std::views::zip(result.buffers_, result.dataset_info_.component_info)) {
+            Idx const size = component_info.elements_per_scenario * batch_size;
+            component_info.total_elements = size;
+        }
+        return result;
+    }
+
   private:
     MetaData const* meta_data_;
     DatasetInfo dataset_info_;
diff --git a/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/basics.h b/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/basics.h
index 1986f35dd7..9153d62d88 100644
--- a/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/basics.h
+++ b/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/basics.h
@@ -137,6 +137,12 @@ typedef struct PGM_WritableDataset PGM_WritableDataset;
  * @brief Opaque struct for the information of the dataset.
  */
 typedef struct PGM_DatasetInfo PGM_DatasetInfo;
+
+/**
+ * @brief Opaque struct for the multi dimensional dataset class.
+ *
+ */
+typedef struct PGM_MultiDimensionalDataset PGM_MultiDimensionalDataset;
 #endif
 
 // NOLINTEND(modernize-use-using)
diff --git a/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/dataset.h b/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/dataset.h
index 6edd6dd0d0..6db842f495 100644
--- a/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/dataset.h
+++ b/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/dataset.h
@@ -309,6 +309,36 @@ PGM_API void PGM_dataset_mutable_add_attribute_buffer(PGM_Handle* handle, PGM_Mu
  */
 PGM_API PGM_DatasetInfo const* PGM_dataset_mutable_get_info(PGM_Handle* handle, PGM_MutableDataset const* dataset);
 
+/**
+ * @brief Create a PGM_MultiDimensionalDataset from multiple PGM_ConstDataset instances
+ *
+ * @param handle
+ * @param const_datasets
+ * @param n_datasets
+ * @return
+ */
+PGM_API PGM_MultiDimensionalDataset*
+PGM_dataset_create_multidimensional_from_const(PGM_Handle* handle, PGM_ConstDataset const** const_datasets,
+                                               PGM_Idx n_datasets);
+
+/**
+ * @brief Get the array pointer from a PGM_MultiDimensionalDataset
+ *
+ * @param handle
+ * @param multidimensional_dataset
+ * @return
+ */
+PGM_API PGM_ConstDataset const*
+PGM_get_array_pointer_from_multidimensional(PGM_Handle* handle,
+                                            PGM_MultiDimensionalDataset const* multidimensional_dataset);
+
+/**
+ * @brief destroy the multidimensional dataset object
+ *
+ * @param multidimensional_dataset
+ */
+PGM_API void PGM_destroy_multidimensional_dataset(PGM_MultiDimensionalDataset* multidimensional_dataset);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/options.h b/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/options.h
index 3671727a00..867a3c517c 100644
--- a/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/options.h
+++ b/power_grid_model_c/power_grid_model_c/include/power_grid_model_c/options.h
@@ -121,6 +121,15 @@ PGM_API void PGM_set_short_circuit_voltage_scaling(PGM_Handle* handle, PGM_Optio
  */
 PGM_API void PGM_set_tap_changing_strategy(PGM_Handle* handle, PGM_Options* opt, PGM_Idx tap_changing_strategy);
 
+/**
+ * @brief Specify the batch dimension for batch calculations
+ *
+ * @param handle
+ * @param opt pointer to option instance
+ * @param batch_dimension dimension of the batch calculation
+ */
+PGM_API void PGM_set_batch_dimension(PGM_Handle* handle, PGM_Options* opt, PGM_Idx batch_dimension);
+
 /**
  * @brief Enable/disable experimental features.
  *
diff --git a/power_grid_model_c/power_grid_model_c/src/dataset.cpp b/power_grid_model_c/power_grid_model_c/src/dataset.cpp
index 0b9e3bee6e..32f785075d 100644
--- a/power_grid_model_c/power_grid_model_c/src/dataset.cpp
+++ b/power_grid_model_c/power_grid_model_c/src/dataset.cpp
@@ -160,3 +160,28 @@ void PGM_dataset_mutable_add_attribute_buffer(PGM_Handle* handle, PGM_MutableDat
 PGM_DatasetInfo const* PGM_dataset_mutable_get_info(PGM_Handle* /*unused*/, PGM_MutableDataset const* dataset) {
     return &dataset->get_description();
 }
+
+PGM_MultiDimensionalDataset* PGM_dataset_create_multidimensional_from_const(PGM_Handle* handle,
+                                                                            PGM_ConstDataset const** const_datasets,
+                                                                            PGM_Idx n_datasets) {
+    return call_with_catch(
+        handle,
+        [const_datasets, n_datasets]() {
+            auto* multidimensional_dataset = new PGM_MultiDimensionalDataset();
+            for (PGM_Idx i = 0; i < n_datasets; ++i) {
+                multidimensional_dataset->emplace_back(*const_datasets[i]);
+            }
+            return multidimensional_dataset;
+        },
+        PGM_regular_error);
+}
+
+PGM_ConstDataset const*
+PGM_get_array_pointer_from_multidimensional(PGM_Handle* /*unused*/,
+                                            PGM_MultiDimensionalDataset const* multidimensional_dataset) {
+    return multidimensional_dataset->data();
+}
+
+void PGM_destroy_multidimensional_dataset(PGM_MultiDimensionalDataset* multidimensional_dataset) {
+    delete multidimensional_dataset;
+}
diff --git a/power_grid_model_c/power_grid_model_c/src/forward_declarations.hpp b/power_grid_model_c/power_grid_model_c/src/forward_declarations.hpp
index f68d7def3f..12532cd5f7 100644
--- a/power_grid_model_c/power_grid_model_c/src/forward_declarations.hpp
+++ b/power_grid_model_c/power_grid_model_c/src/forward_declarations.hpp
@@ -10,6 +10,8 @@
 
 #include <power_grid_model/auxiliary/dataset_fwd.hpp>
 
+#include <vector>
+
 // forward declare all referenced struct/class in C++ core
 // alias them in the root namespace
 
@@ -36,3 +38,4 @@ using PGM_ConstDataset = power_grid_model::meta_data::Dataset<power_grid_model::
 using PGM_MutableDataset = power_grid_model::meta_data::Dataset<power_grid_model::mutable_dataset_t>;
 using PGM_WritableDataset = power_grid_model::meta_data::Dataset<power_grid_model::writable_dataset_t>;
 using PGM_DatasetInfo = power_grid_model::meta_data::DatasetInfo;
+using PGM_MultiDimensionalDataset = std::vector<PGM_ConstDataset>;
diff --git a/power_grid_model_c/power_grid_model_c/src/model.cpp b/power_grid_model_c/power_grid_model_c/src/model.cpp
index 8310c62eae..840e3cf64d 100644
--- a/power_grid_model_c/power_grid_model_c/src/model.cpp
+++ b/power_grid_model_c/power_grid_model_c/src/model.cpp
@@ -15,6 +15,9 @@
 #include <power_grid_model/common/common.hpp>
 #include <power_grid_model/main_model.hpp>
 
+#include <memory>
+#include <numeric>
+
 namespace {
 using namespace power_grid_model;
 } // namespace
@@ -55,6 +58,7 @@ void PGM_get_indexer(PGM_Handle* handle, PGM_PowerGridModel const* model, char c
         PGM_regular_error);
 }
 
+// helper functions
 namespace {
 void check_no_experimental_features_used(MainModel const& model, MainModel::Options const& opt) {
     // optionally add experimental feature checks here
@@ -142,9 +146,11 @@ constexpr auto extract_calculation_options(PGM_Options const& opt) {
 }
 } // namespace
 
-// run calculation
-void PGM_calculate(PGM_Handle* handle, PGM_PowerGridModel* model, PGM_Options const* opt,
-                   PGM_MutableDataset const* output_dataset, PGM_ConstDataset const* batch_dataset) {
+// calculation implementation
+namespace {
+
+void calculate_impl(PGM_Handle* handle, PGM_PowerGridModel* model, PGM_Options const* opt,
+                    PGM_MutableDataset const* output_dataset, PGM_ConstDataset const* batch_dataset) {
     PGM_clear_error(handle);
     // check dataset integrity
     if ((batch_dataset != nullptr) && (!batch_dataset->is_batch() || !output_dataset->is_batch())) {
@@ -180,5 +186,83 @@ void PGM_calculate(PGM_Handle* handle, PGM_PowerGridModel* model, PGM_Options co
     }
 }
 
+void merge_batch_error_msgs(PGM_Handle* handle, PGM_Handle const& local_handle, Idx scenario_offset, Idx stride_size) {
+    if (local_handle.err_code == PGM_no_error) {
+        return;
+    }
+    handle->err_code = PGM_batch_error;
+    if (local_handle.err_code == PGM_batch_error) {
+        for (auto&& [idx, err_msg] : std::views::zip(local_handle.failed_scenarios, local_handle.batch_errs)) {
+            handle->failed_scenarios.push_back(idx + scenario_offset);
+            handle->batch_errs.push_back(err_msg);
+        }
+    } else {
+        for (Idx i = 0; i < stride_size; ++i) {
+            handle->failed_scenarios.push_back(scenario_offset + i);
+            handle->batch_errs.push_back(local_handle.err_msg);
+        }
+    }
+}
+
+} // namespace
+
+// run calculation
+void PGM_calculate(PGM_Handle* handle, PGM_PowerGridModel* model, PGM_Options const* opt,
+                   PGM_MutableDataset const* output_dataset, PGM_ConstDataset const* batch_dataset) {
+    // if dimension is zero, no batch calculation, force pointer to NULL
+    if (opt->batch_dimension == 0) {
+        batch_dataset = nullptr;
+    }
+
+    // for dimensions which are 1D batch or default (-1), call implementation directly
+    if (opt->batch_dimension < 2) {
+        calculate_impl(handle, model, opt, output_dataset, batch_dataset);
+        return;
+    }
+
+    // get stride size of the rest of dimensions
+    Idx const first_batch_size = batch_dataset->batch_size();
+    Idx const stride_size =
+        std::transform_reduce(batch_dataset + 1, batch_dataset + opt->batch_dimension, Idx{1}, std::multiplies{},
+                              [](PGM_ConstDataset const& ds) { return ds.batch_size(); });
+
+    // loop over the first dimension batche
+    for (Idx i = 0; i < first_batch_size; ++i) {
+        // a new handle
+        PGM_Handle local_handle{};
+        // deep opt is one dimension less
+        PGM_Options deep_opt = *opt;
+        --deep_opt.batch_dimension;
+        // create sliced datasets for the rest of dimensions
+        PGM_ConstDataset const single_update_dataset = batch_dataset->get_individual_scenario(i);
+        PGM_MutableDataset const sliced_output_dataset =
+            output_dataset->get_slice_scenario(i * stride_size, (i + 1) * stride_size);
+
+        // create a model copy
+        std::unique_ptr<PGM_PowerGridModel> const local_model{PGM_copy_model(&local_handle, model)};
+        if (local_handle.err_code != PGM_no_error) {
+            merge_batch_error_msgs(handle, local_handle, i * stride_size, stride_size);
+            continue;
+        }
+
+        // apply the update
+        PGM_update_model(&local_handle, local_model.get(), &single_update_dataset);
+        if (local_handle.err_code != PGM_no_error) {
+            merge_batch_error_msgs(handle, local_handle, i * stride_size, stride_size);
+            continue;
+        }
+
+        // if the deep opt have less than 2 dimensions, call implementation directly
+        if (deep_opt.batch_dimension < 2) {
+            calculate_impl(&local_handle, local_model.get(), &deep_opt, &sliced_output_dataset, batch_dataset + 1);
+        } else {
+            // recursive call
+            PGM_calculate(&local_handle, local_model.get(), &deep_opt, &sliced_output_dataset, batch_dataset + 1);
+        }
+        // merge errors
+        merge_batch_error_msgs(handle, local_handle, i * stride_size, stride_size);
+    }
+}
+
 // destroy model
 void PGM_destroy_model(PGM_PowerGridModel* model) { delete model; }
diff --git a/power_grid_model_c/power_grid_model_c/src/options.cpp b/power_grid_model_c/power_grid_model_c/src/options.cpp
index 5c22d7acaa..15b517f757 100644
--- a/power_grid_model_c/power_grid_model_c/src/options.cpp
+++ b/power_grid_model_c/power_grid_model_c/src/options.cpp
@@ -33,6 +33,9 @@ void PGM_set_short_circuit_voltage_scaling(PGM_Handle* /* handle */, PGM_Options
 void PGM_set_tap_changing_strategy(PGM_Handle* /* handle */, PGM_Options* opt, PGM_Idx tap_changing_strategy) {
     opt->tap_changing_strategy = tap_changing_strategy;
 }
+void PGM_set_batch_dimension(PGM_Handle* /* handle */, PGM_Options* opt, PGM_Idx batch_dimension) {
+    opt->batch_dimension = batch_dimension;
+}
 void PGM_set_experimental_features(PGM_Handle* /* handle */, PGM_Options* opt, PGM_Idx experimental_features) {
     opt->experimental_features = experimental_features;
 }
diff --git a/power_grid_model_c/power_grid_model_c/src/options.hpp b/power_grid_model_c/power_grid_model_c/src/options.hpp
index 37da22cb77..b72b4d2d5f 100644
--- a/power_grid_model_c/power_grid_model_c/src/options.hpp
+++ b/power_grid_model_c/power_grid_model_c/src/options.hpp
@@ -25,5 +25,6 @@ struct PGM_Options {
     Idx threading{-1};
     Idx short_circuit_voltage_scaling{PGM_short_circuit_voltage_scaling_maximum};
     Idx tap_changing_strategy{PGM_tap_changing_strategy_disabled};
+    Idx batch_dimension{-1};
     Idx experimental_features{PGM_experimental_features_disabled};
 };
diff --git a/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/model.hpp b/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/model.hpp
index 822e383ae0..d2e3cd38cc 100644
--- a/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/model.hpp
+++ b/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/model.hpp
@@ -13,6 +13,8 @@
 
 #include "power_grid_model_c/model.h"
 
+#include <vector>
+
 namespace power_grid_model_cpp {
 class Model {
   public:
@@ -54,6 +56,25 @@ class Model {
         handle_.call_with(PGM_calculate, get(), opt.get(), output_dataset.get(), nullptr);
     }
 
+    void calculate(Options const& opt, DatasetMutable const& output_dataset,
+                   std::vector<DatasetConst> const& batch_datasets) {
+        // create multidimensional dataset from the span of datasets
+        std::vector<PGM_ConstDataset const*> dataset_ptrs;
+        dataset_ptrs.reserve(batch_datasets.size());
+        for (auto const& ds : batch_datasets) {
+            dataset_ptrs.push_back(ds.get());
+        }
+        auto multidimensional_dataset =
+            detail::UniquePtr<PGM_MultiDimensionalDataset, &PGM_destroy_multidimensional_dataset>{
+                handle_.call_with(PGM_dataset_create_multidimensional_from_const, dataset_ptrs.data(),
+                                  static_cast<Idx>(dataset_ptrs.size()))};
+        RawConstDataset const* batch_dataset_array_pointer =
+            PGM_get_array_pointer_from_multidimensional(nullptr, multidimensional_dataset.get());
+
+        // call calculate with the multidimensional dataset
+        handle_.call_with(PGM_calculate, get(), opt.get(), output_dataset.get(), batch_dataset_array_pointer);
+    }
+
   private:
     Handle handle_{};
     detail::UniquePtr<PowerGridModel, &PGM_destroy_model> model_;
diff --git a/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/options.hpp b/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/options.hpp
index 5603650e5f..3b29364546 100644
--- a/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/options.hpp
+++ b/power_grid_model_c/power_grid_model_cpp/include/power_grid_model_cpp/options.hpp
@@ -39,6 +39,10 @@ class Options {
         handle_.call_with(PGM_set_tap_changing_strategy, get(), tap_changing_strategy);
     }
 
+    void set_batch_dimension(Idx batch_dimension) {
+        handle_.call_with(PGM_set_batch_dimension, get(), batch_dimension);
+    }
+
     void set_experimental_features(Idx experimental_features) {
         handle_.call_with(PGM_set_experimental_features, get(), experimental_features);
     }
diff --git a/src/power_grid_model/_core/options.py b/src/power_grid_model/_core/options.py
index f8f5a910ba..1bb9aa419e 100644
--- a/src/power_grid_model/_core/options.py
+++ b/src/power_grid_model/_core/options.py
@@ -45,6 +45,7 @@ class Options:
     tap_changing_strategy = OptionSetter(pgc.set_tap_changing_strategy)
     short_circuit_voltage_scaling = OptionSetter(pgc.set_short_circuit_voltage_scaling)
     experimental_features = OptionSetter(pgc.set_experimental_features)
+    batch_dimension = OptionSetter(pgc.set_batch_dimension)
 
     @property
     def opt(self) -> OptionsPtr:
diff --git a/src/power_grid_model/_core/power_grid_core.py b/src/power_grid_model/_core/power_grid_core.py
index c276303273..3a03a34d49 100644
--- a/src/power_grid_model/_core/power_grid_core.py
+++ b/src/power_grid_model/_core/power_grid_core.py
@@ -87,6 +87,12 @@ class WritableDatasetPtr(c_void_p):
     """
 
 
+class MultiDimensionDatasetPtr(c_void_p):
+    """
+    Pointer to multi-dimensional dataset
+    """
+
+
 class DatasetInfoPtr(c_void_p):
     """
     Pointer to dataset info
@@ -338,6 +344,10 @@ def set_max_iter(self, opt: OptionsPtr, max_iter: int) -> None:  # type: ignore[
     def set_threading(self, opt: OptionsPtr, threading: int) -> None:  # type: ignore[empty-body]
         pass  # pragma: no cover
 
+    @make_c_binding
+    def set_batch_dimension(self, opt: OptionsPtr, batch_dimension: int) -> None:  # type: ignore[empty-body]
+        pass  # pragma: no cover
+
     @make_c_binding
     def create_model(  # type: ignore[empty-body]
         self,
@@ -505,6 +515,28 @@ def dataset_writable_set_attribute_buffer(
     ) -> None:  # type: ignore[empty-body]
         pass  # pragma: no cover
 
+    @make_c_binding
+    def dataset_create_multidimensional_from_const(  # type: ignore[empty-body]
+        self,
+        const_datasets: POINTER(ConstDatasetPtr),  # type: ignore[valid-type]
+        n_datasets: int,
+    ) -> MultiDimensionDatasetPtr:
+        pass  # pragma: no cover
+
+    @make_c_binding
+    def get_array_pointer_from_multidimensional(  # type: ignore[empty-body]
+        self,
+        multidimensional_dataset: MultiDimensionDatasetPtr,
+    ) -> ConstDatasetPtr:  # type: ignore[valid-type]
+        pass  # pragma: no cover
+
+    @make_c_binding
+    def destroy_multidimensional_dataset(  # type: ignore[empty-body]
+        self,
+        multidimensional_dataset: MultiDimensionDatasetPtr,
+    ) -> None:
+        pass  # pragma: no cover
+
     @make_c_binding
     def create_deserializer_from_binary_buffer(  # type: ignore[empty-body]
         self, data: bytes, size: int, serialization_format: int
diff --git a/src/power_grid_model/_core/power_grid_dataset.py b/src/power_grid_model/_core/power_grid_dataset.py
index f6013f16f0..88bbf6c359 100644
--- a/src/power_grid_model/_core/power_grid_dataset.py
+++ b/src/power_grid_model/_core/power_grid_dataset.py
@@ -32,6 +32,7 @@
 from power_grid_model._core.power_grid_core import (
     ConstDatasetPtr,
     DatasetInfoPtr,
+    MultiDimensionDatasetPtr,
     MutableDatasetPtr,
     WritableDatasetPtr,
     power_grid_core as pgc,
@@ -354,6 +355,57 @@ def __del__(self):
         pgc.destroy_dataset_const(self._const_dataset)
 
 
+class CMultiDimensionalDataset:
+    """
+    A proxy of a multi-dimensional dataset representing multiple const dataset for multidimensional batch scenarios.
+    """
+
+    _multi_dimensional_dataset: MultiDimensionDatasetPtr
+    _const_datasets: list[CConstDataset]
+    _total_batch_size: int
+
+    def __new__(cls, datasets: list[CConstDataset]):
+        instance = super().__new__(cls)
+        instance._multi_dimensional_dataset = MultiDimensionDatasetPtr()
+        instance._const_datasets = datasets
+
+        DatasetsPtrArray = ConstDatasetPtr * len(datasets)
+        datasets_ptr_array = DatasetsPtrArray()
+        total_batch_size = 1
+        for idx, dataset in enumerate(datasets):
+            datasets_ptr_array[idx] = dataset.get_dataset_ptr()
+            total_batch_size *= dataset.get_info().batch_size()
+        instance._total_batch_size = total_batch_size
+
+        instance._multi_dimensional_dataset = pgc.dataset_create_multidimensional_from_const(
+            datasets_ptr_array, len(datasets)
+        )
+        assert_no_error()
+
+        return instance
+
+    def get_array_ptr(self) -> ConstDatasetPtr:
+        """
+        Get the raw underlying multi-dimensional dataset pointer.
+
+        Returns:
+            MultiDimensionDatasetPtr: the raw underlying multi-dimensional dataset pointer.
+        """
+        return pgc.get_array_pointer_from_multidimensional(self._multi_dimensional_dataset)
+
+    def get_total_batch_size(self) -> int:
+        """
+        Get the total batch size of the multi-dimensional dataset.
+
+        Returns:
+            int: the total batch size of the multi-dimensional dataset.
+        """
+        return self._total_batch_size
+
+    def __del__(self):
+        pgc.destroy_multidimensional_dataset(self._multi_dimensional_dataset)
+
+
 class CWritableDataset:
     """
     A view of a Power Grid Model-owned dataset.
diff --git a/src/power_grid_model/_core/power_grid_model.py b/src/power_grid_model/_core/power_grid_model.py
index f182d04c4c..2f9c0d2803 100644
--- a/src/power_grid_model/_core/power_grid_model.py
+++ b/src/power_grid_model/_core/power_grid_model.py
@@ -49,6 +49,7 @@
 from power_grid_model._core.index_integer import IdNp, IdxNp
 from power_grid_model._core.options import Options
 from power_grid_model._core.power_grid_core import ConstDatasetPtr, IDPtr, IdxPtr, ModelPtr, power_grid_core as pgc
+from power_grid_model._core.power_grid_dataset import CMultiDimensionalDataset
 from power_grid_model._core.typing import ComponentAttributeMapping, ComponentAttributeMappingDict
 
 
@@ -261,7 +262,7 @@ def _calculate_impl(  # noqa: PLR0913
         self,
         calculation_type: CalculationType,
         symmetric: bool,
-        update_data: Dataset | None,
+        update_data: Dataset | list[Dataset] | None,
         output_component_types: ComponentAttributeMapping,
         options: Options,
         continue_on_batch_error: bool,
@@ -283,15 +284,19 @@ def _calculate_impl(  # noqa: PLR0913
         Returns:
         """
         self._batch_error = None
-        is_batch = update_data is not None
-
-        if update_data is not None:
-            prepared_update = prepare_update_view(update_data)
-            update_ptr = prepared_update.get_dataset_ptr()
-            batch_size = prepared_update.get_info().batch_size()
+        if update_data is None:
+            is_batch = False
+            update_data = []
         else:
-            update_ptr = ConstDatasetPtr()
-            batch_size = 1
+            is_batch = True
+        if not isinstance(update_data, list):
+            update_data = [update_data]
+        options.batch_dimension = len(update_data)
+        update_data = [_map_to_component_types(x) for x in update_data]
+
+        prepared_update = CMultiDimensionalDataset([prepare_update_view(x) for x in update_data])
+        update_ptr: ConstDatasetPtr = prepared_update.get_array_ptr()
+        batch_size = prepared_update.get_total_batch_size()
 
         output_data = self._construct_output(
             output_component_types=output_component_types,
@@ -329,7 +334,7 @@ def _calculate_power_flow(  # noqa: PLR0913
         error_tolerance: float = 1e-8,
         max_iterations: int = 20,
         calculation_method: CalculationMethod | str = CalculationMethod.newton_raphson,
-        update_data: Dataset | None = None,
+        update_data: Dataset | list[Dataset] | None = None,
         threading: int = -1,
         output_component_types: ComponentAttributeMapping = None,
         continue_on_batch_error: bool = False,
@@ -366,7 +371,7 @@ def _calculate_state_estimation(  # noqa: PLR0913
         error_tolerance: float = 1e-8,
         max_iterations: int = 20,
         calculation_method: CalculationMethod | str = CalculationMethod.iterative_linear,
-        update_data: Dataset | None = None,
+        update_data: Dataset | list[Dataset] | None = None,
         threading: int = -1,
         output_component_types: ComponentAttributeMapping = None,
         continue_on_batch_error: bool = False,
@@ -398,7 +403,7 @@ def _calculate_short_circuit(  # noqa: PLR0913
         self,
         *,
         calculation_method: CalculationMethod | str = CalculationMethod.iec60909,
-        update_data: Dataset | None = None,
+        update_data: Dataset | list[Dataset] | None = None,
         threading: int = -1,
         output_component_types: ComponentAttributeMapping = None,
         continue_on_batch_error: bool = False,
@@ -495,7 +500,7 @@ def calculate_power_flow(
         error_tolerance: float = ...,
         max_iterations: int = ...,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: None | set[ComponentTypeVar] | list[ComponentTypeVar] = ...,
         continue_on_batch_error: bool = ...,
@@ -510,7 +515,7 @@ def calculate_power_flow(
         error_tolerance: float = ...,
         max_iterations: int = ...,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: ComponentAttributeFilterOptions = ...,
         continue_on_batch_error: bool = ...,
@@ -525,7 +530,7 @@ def calculate_power_flow(
         error_tolerance: float = ...,
         max_iterations: int = ...,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: ComponentAttributeMappingDict = ...,
         continue_on_batch_error: bool = ...,
@@ -539,7 +544,7 @@ def calculate_power_flow(  # noqa: PLR0913
         error_tolerance: float = 1e-8,
         max_iterations: int = 20,
         calculation_method: CalculationMethod | str = CalculationMethod.newton_raphson,
-        update_data: BatchDataset | None = None,
+        update_data: BatchDataset | list[BatchDataset] | None = None,
         threading: int = -1,
         output_component_types: ComponentAttributeMapping = None,
         continue_on_batch_error: bool = False,
@@ -623,7 +628,7 @@ def calculate_power_flow(  # noqa: PLR0913
             error_tolerance=error_tolerance,
             max_iterations=max_iterations,
             calculation_method=calculation_method,
-            update_data=(_map_to_component_types(update_data) if update_data is not None else None),
+            update_data=update_data,
             threading=threading,
             output_component_types=output_component_types,
             continue_on_batch_error=continue_on_batch_error,
@@ -681,7 +686,7 @@ def calculate_state_estimation(
         error_tolerance: float = ...,
         max_iterations: int = ...,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: None | set[ComponentTypeVar] | list[ComponentTypeVar] = ...,
         continue_on_batch_error: bool = ...,
@@ -695,7 +700,7 @@ def calculate_state_estimation(
         error_tolerance: float = ...,
         max_iterations: int = ...,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: ComponentAttributeFilterOptions = ...,
         continue_on_batch_error: bool = ...,
@@ -709,7 +714,7 @@ def calculate_state_estimation(
         error_tolerance: float = ...,
         max_iterations: int = ...,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: ComponentAttributeMappingDict = ...,
         continue_on_batch_error: bool = ...,
@@ -722,7 +727,7 @@ def calculate_state_estimation(  # noqa: PLR0913
         error_tolerance: float = 1e-8,
         max_iterations: int = 20,
         calculation_method: CalculationMethod | str = CalculationMethod.iterative_linear,
-        update_data: BatchDataset | None = None,
+        update_data: BatchDataset | list[BatchDataset] | None = None,
         threading: int = -1,
         output_component_types: ComponentAttributeMapping = None,
         continue_on_batch_error: bool = False,
@@ -802,7 +807,7 @@ def calculate_state_estimation(  # noqa: PLR0913
             error_tolerance=error_tolerance,
             max_iterations=max_iterations,
             calculation_method=calculation_method,
-            update_data=(_map_to_component_types(update_data) if update_data is not None else None),
+            update_data=update_data,
             threading=threading,
             output_component_types=output_component_types,
             continue_on_batch_error=continue_on_batch_error,
@@ -850,7 +855,7 @@ def calculate_short_circuit(
         self,
         *,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: None | set[ComponentTypeVar] | list[ComponentTypeVar] = ...,
         continue_on_batch_error: bool = ...,
@@ -862,7 +867,7 @@ def calculate_short_circuit(
         self,
         *,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: ComponentAttributeFilterOptions = ...,
         continue_on_batch_error: bool = ...,
@@ -874,7 +879,7 @@ def calculate_short_circuit(
         self,
         *,
         calculation_method: CalculationMethod | str = ...,
-        update_data: BatchDataset = ...,
+        update_data: BatchDataset | list[BatchDataset] = ...,
         threading: int = ...,
         output_component_types: ComponentAttributeMappingDict = ...,
         continue_on_batch_error: bool = ...,
@@ -885,7 +890,7 @@ def calculate_short_circuit(  # noqa: PLR0913
         self,
         *,
         calculation_method: CalculationMethod | str = CalculationMethod.iec60909,
-        update_data: BatchDataset | None = None,
+        update_data: BatchDataset | list[BatchDataset] | None = None,
         threading: int = -1,
         output_component_types: ComponentAttributeMapping = None,
         continue_on_batch_error: bool = False,
@@ -957,7 +962,7 @@ def calculate_short_circuit(  # noqa: PLR0913
         """
         return self._calculate_short_circuit(
             calculation_method=calculation_method,
-            update_data=(_map_to_component_types(update_data) if update_data is not None else None),
+            update_data=update_data,
             threading=threading,
             output_component_types=output_component_types,
             continue_on_batch_error=continue_on_batch_error,
diff --git a/tests/native_api_tests/CMakeLists.txt b/tests/native_api_tests/CMakeLists.txt
index e58c80b726..739b8be3cc 100644
--- a/tests/native_api_tests/CMakeLists.txt
+++ b/tests/native_api_tests/CMakeLists.txt
@@ -10,6 +10,7 @@ set(PROJECT_SOURCES
     "test_api_model_update.cpp"
     "test_api_serialization.cpp"
     "test_api_utils.cpp"
+    "test_api_model_multi_dimension.cpp"
 )
 
 add_executable(power_grid_model_api_tests ${PROJECT_SOURCES})
diff --git a/tests/native_api_tests/test_api_model_multi_dimension.cpp b/tests/native_api_tests/test_api_model_multi_dimension.cpp
new file mode 100644
index 0000000000..b9370426ae
--- /dev/null
+++ b/tests/native_api_tests/test_api_model_multi_dimension.cpp
@@ -0,0 +1,109 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+#include "load_dataset.hpp"
+
+#include "power_grid_model_cpp.hpp"
+
+#include <power_grid_model_c/dataset_definitions.h>
+
+#include <doctest/doctest.h>
+
+#include <array>
+#include <complex>
+#include <numbers>
+#include <string>
+#include <vector>
+
+namespace power_grid_model_cpp {
+namespace {
+using namespace std::string_literals;
+using power_grid_model_cpp_test::load_dataset;
+using std::numbers::sqrt3;
+
+// input
+auto const complete_state_json = R"json({
+  "version": "1.0",
+  "type": "input",
+  "is_batch": false,
+  "attributes": {},
+  "data": {
+    "sym_load": [
+      {"id": 2, "node": 0, "status": 1, "type": 0, "p_specified": 0, "q_specified": 0}
+    ],
+    "source": [
+      {"id": 1, "node": 0, "status": 1, "u_ref": 1, "sk": 1e20}
+    ],
+    "node": [
+      {"id": 0, "u_rated": 10e3}
+    ]
+  }
+})json"s;
+
+} // namespace
+
+TEST_CASE("API Model Multi-Dimension") {
+    // model
+    auto const owning_input_dataset = load_dataset(complete_state_json);
+    auto const& input_dataset = owning_input_dataset.dataset;
+    Model model{50.0, input_dataset};
+
+    // 3-D batch update
+    double const u_rated = 10e3;
+    std::vector<double> const u_ref{0.9, 1.0, 1.1};
+    std::vector<double> const p_specified{1e6, 2e6, 3e6, 4e6};
+    std::vector<double> const q_specified{0.1e6, 0.2e6, 0.3e6, 0.4e6, 0.5e6};
+    Idx const size_u_ref = static_cast<Idx>(u_ref.size());
+    Idx const size_p_specified = static_cast<Idx>(p_specified.size());
+    Idx const size_q_specified = static_cast<Idx>(q_specified.size());
+    Idx const total_batch_size = size_u_ref * size_p_specified * size_q_specified;
+
+    // calculate source current manually
+    std::vector<double> i_source_ref(total_batch_size);
+    for (Idx i = 0; i < size_u_ref; ++i) {
+        for (Idx j = 0; j < size_p_specified; ++j) {
+            for (Idx k = 0; k < size_q_specified; ++k) {
+                Idx const index = i * size_p_specified * size_q_specified + j * size_q_specified + k;
+                double const s = std::abs(std::complex<double>{p_specified[j], q_specified[k]});
+                i_source_ref[index] = s / (sqrt3 * u_rated * u_ref[i]);
+            }
+        }
+    }
+
+    // construct batch update dataset
+    std::vector<DatasetConst> batch_datasets;
+    batch_datasets.emplace_back("update", true, size_u_ref);
+    batch_datasets.emplace_back("update", true, size_p_specified);
+    batch_datasets.emplace_back("update", true, size_q_specified);
+
+    DatasetConst& batch_u_ref = batch_datasets[0];
+    batch_u_ref.add_buffer("source", 1, size_u_ref, nullptr, nullptr);
+    batch_u_ref.add_attribute_buffer("source", "u_ref", u_ref.data());
+    DatasetConst& batch_p_specified = batch_datasets[1];
+    batch_p_specified.add_buffer("sym_load", 1, size_p_specified, nullptr, nullptr);
+    batch_p_specified.add_attribute_buffer("sym_load", "p_specified", p_specified.data());
+    DatasetConst& batch_q_specified = batch_datasets[2];
+    batch_q_specified.add_buffer("sym_load", 1, size_q_specified, nullptr, nullptr);
+    batch_q_specified.add_attribute_buffer("sym_load", "q_specified", q_specified.data());
+
+    // output dataset
+    std::vector<double> i_source_result(total_batch_size);
+    DatasetMutable batch_output_dataset{"sym_output", true, total_batch_size};
+    batch_output_dataset.add_buffer("source", 1, total_batch_size, nullptr, nullptr);
+    batch_output_dataset.add_attribute_buffer("source", "i", i_source_result.data());
+
+    // options
+    Options options{};
+    options.set_batch_dimension(3);
+
+    // calculate
+    model.calculate(options, batch_output_dataset, batch_datasets);
+
+    // check results
+    for (Idx idx = 0; idx < total_batch_size; ++idx) {
+        CHECK(i_source_result[idx] == doctest::Approx(i_source_ref[idx]));
+    }
+}
+
+} // namespace power_grid_model_cpp
diff --git a/tests/unit/test_multi_dimensional_batch.py b/tests/unit/test_multi_dimensional_batch.py
new file mode 100644
index 0000000000..2ca9e13faa
--- /dev/null
+++ b/tests/unit/test_multi_dimensional_batch.py
@@ -0,0 +1,48 @@
+# SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+#
+# SPDX-License-Identifier: MPL-2.0
+
+import json
+
+import numpy as np
+
+from power_grid_model import PowerGridModel
+from power_grid_model.utils import json_deserialize
+
+input_data = {
+    "version": "1.0",
+    "type": "input",
+    "is_batch": False,
+    "attributes": {},
+    "data": {
+        "sym_load": [{"id": 2, "node": 0, "status": 1, "type": 0, "p_specified": 0, "q_specified": 0}],
+        "source": [{"id": 1, "node": 0, "status": 1, "u_ref": 1, "sk": 1e20}],
+        "node": [{"id": 0, "u_rated": 10e3}],
+    },
+}
+
+input_data_json = json.dumps(input_data)
+
+
+def test_multi_dimensional_batch():
+    input_dataset = json_deserialize(input_data_json)
+    pgm = PowerGridModel(input_dataset)
+
+    u_rated = 10e3
+    u_ref = np.array([0.9, 1.0, 1.1], dtype=np.float64).reshape(-1, 1)
+    p_specified = np.array([1e6, 2e6, 3e6, 4e6], dtype=np.float64).reshape(-1, 1)
+    q_specified = np.array([0.1e6, 0.2e6, 0.3e6, 0.4e6, 0.5e6], dtype=np.float64).reshape(-1, 1)
+    i_source_ref = np.abs(p_specified.reshape(1, -1, 1) + 1j * q_specified.reshape(1, 1, -1)) / (
+        u_ref.reshape(-1, 1, 1) * u_rated * np.sqrt(3)
+    )
+    i_source_ref = i_source_ref.ravel()
+
+    u_ref_batch = {"source": {"u_ref": u_ref}}
+    p_specified_batch = {"sym_load": {"p_specified": p_specified}}
+    q_specified_batch = {"sym_load": {"q_specified": q_specified}}
+
+    result = pgm.calculate_power_flow(
+        update_data=[u_ref_batch, p_specified_batch, q_specified_batch], output_component_types={"source": ["i"]}
+    )
+
+    assert np.allclose(result["source"]["i"].ravel(), i_source_ref)