diff --git a/src/language/PEGGrammar.hpp b/src/language/PEGGrammar.hpp
index 20ec81befa009efd926820b0a6333be0cc48dc5e..beae437f0e5d0f9231b51afc0671e29e7ecfef30 100644
--- a/src/language/PEGGrammar.hpp
+++ b/src/language/PEGGrammar.hpp
@@ -181,7 +181,7 @@ struct postfix_operator : seq< sor< post_plusplus, post_minusminus>, ignored > {
 struct open_bracket : seq< one< '[' >, ignored > {};
 struct close_bracket : seq< one< ']' >, ignored > {};
 
-struct subscript_expression : if_must< open_bracket, expression, close_bracket >{};
+struct subscript_expression : if_must< open_bracket, list_must<expression, COMMA>, close_bracket >{};
 
 struct postfix_expression : seq< primary_expression, star< sor< subscript_expression , postfix_operator> > >{};
 
diff --git a/src/language/ast/ASTBuilder.cpp b/src/language/ast/ASTBuilder.cpp
index f9e3edceb480b01c608997e31c82fb133b8e7c35..ccf0f42db5c4144ef91475d3fa9c3c24520d1e35 100644
--- a/src/language/ast/ASTBuilder.cpp
+++ b/src/language/ast/ASTBuilder.cpp
@@ -108,21 +108,19 @@ struct ASTBuilder::simplify_unary : parse_tree::apply<ASTBuilder::simplify_unary
     }
 
     if (n->is_type<language::unary_expression>() or n->is_type<language::name_subscript_expression>()) {
-      const size_t child_nb = n->children.size();
-      if (child_nb > 1) {
+      if (n->children.size() > 1) {
         if (n->children[1]->is_type<language::subscript_expression>()) {
-          auto expression = std::move(n->children[0]);
-          for (size_t i = 0; i < child_nb - 1; ++i) {
-            n->children[i] = std::move(n->children[i + 1]);
-          }
+          std::swap(n->children[0], n->children[1]);
 
-          auto& array_subscript_expression = n->children[0];
+          n->children[0]->emplace_back(std::move(n->children[1]));
           n->children.pop_back();
-          array_subscript_expression->children.emplace_back(std::move(expression));
 
-          std::swap(array_subscript_expression->children[0], array_subscript_expression->children[1]);
-
-          array_subscript_expression->m_begin = array_subscript_expression->children[0]->m_begin;
+          auto& array_subscript_expression = n->children[0];
+          const size_t child_nb            = array_subscript_expression->children.size();
+          for (size_t i = 1; i < array_subscript_expression->children.size(); ++i) {
+            std::swap(array_subscript_expression->children[child_nb - i],
+                      array_subscript_expression->children[child_nb - i - 1]);
+          }
 
           transform(n, st...);
         }
diff --git a/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp b/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp
index 350a4070e9a6bc65906627d54b2e2ed87c403cc6..e717044d1dea5fe6b85aec6c41941699c1e5ef65 100644
--- a/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp
@@ -1,5 +1,6 @@
 #include <language/ast/ASTNodeArraySubscriptExpressionBuilder.hpp>
 
+#include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
 #include <language/node_processor/ArraySubscriptProcessor.hpp>
 #include <language/utils/ParseError.hpp>
@@ -27,6 +28,27 @@ ASTNodeArraySubscriptExpressionBuilder::ASTNodeArraySubscriptExpressionBuilder(A
       break;
     }
     }
+  } else if (array_expression.m_data_type == ASTNodeDataType::matrix_t) {
+    Assert(array_expression.m_data_type.nbRows() == array_expression.m_data_type.nbColumns());
+
+    switch (array_expression.m_data_type.nbRows()) {
+    case 1: {
+      node.m_node_processor = std::make_unique<ArraySubscriptProcessor<TinyMatrix<1>>>(node);
+      break;
+    }
+    case 2: {
+      node.m_node_processor = std::make_unique<ArraySubscriptProcessor<TinyMatrix<2>>>(node);
+      break;
+    }
+    case 3: {
+      node.m_node_processor = std::make_unique<ArraySubscriptProcessor<TinyMatrix<3>>>(node);
+      break;
+    }
+    default: {
+      throw ParseError("unexpected error: invalid array dimension", array_expression.begin());
+      break;
+    }
+    }
   } else {
     throw ParseError("unexpected error: invalid array type", array_expression.begin());
   }
diff --git a/src/language/ast/ASTNodeDataTypeBuilder.cpp b/src/language/ast/ASTNodeDataTypeBuilder.cpp
index 1cac7a59f967ea930bcd661c268044a09f778742..2f1612cb2c5becf7d0c9a8755d4627823e3b380e 100644
--- a/src/language/ast/ASTNodeDataTypeBuilder.cpp
+++ b/src/language/ast/ASTNodeDataTypeBuilder.cpp
@@ -523,20 +523,37 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
         throw ParseError(message.str(), n.begin());
       }
     } else if (n.is_type<language::subscript_expression>()) {
-      Assert(n.children.size() == 2, "invalid number of sub-expressions in array subscript expression");
       auto& array_expression = *n.children[0];
-      auto& index_expression = *n.children[1];
 
-      ASTNodeNaturalConversionChecker{index_expression, ASTNodeDataType::build<ASTNodeDataType::int_t>()};
-      if (array_expression.m_data_type != ASTNodeDataType::vector_t) {
+      if (array_expression.m_data_type == ASTNodeDataType::vector_t) {
+        auto& index_expression = *n.children[1];
+        ASTNodeNaturalConversionChecker{index_expression, ASTNodeDataType::build<ASTNodeDataType::int_t>()};
+        n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::double_t>();
+        if (n.children.size() != 2) {
+          std::ostringstream message;
+          message << "invalid index type: " << rang::fgB::yellow << dataTypeName(array_expression.m_data_type)
+                  << rang::style::reset << " requires a single integer";
+          throw ParseError(message.str(), index_expression.begin());
+        }
+      } else if (array_expression.m_data_type == ASTNodeDataType::matrix_t) {
+        for (size_t i = 1; i < n.children.size(); ++i) {
+          ASTNodeNaturalConversionChecker{*n.children[i], ASTNodeDataType::build<ASTNodeDataType::int_t>()};
+        }
+        n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::double_t>();
+
+        if (n.children.size() != 3) {
+          std::ostringstream message;
+          message << "invalid index type: " << rang::fgB::yellow << dataTypeName(n.children[0]->m_data_type)
+                  << rang::style::reset << " requires two integers";
+          throw ParseError(message.str(), n.children[1]->begin());
+        }
+
+      } else {
         std::ostringstream message;
-        message << "invalid types '" << rang::fgB::yellow << dataTypeName(array_expression.m_data_type)
-                << rang::style::reset << '[' << dataTypeName(index_expression.m_data_type) << ']'
-                << "' for array subscript";
+        message << "invalid subscript expression: " << rang::fgB::yellow << dataTypeName(array_expression.m_data_type)
+                << rang::style::reset << " cannot be indexed";
 
         throw ParseError(message.str(), n.begin());
-      } else {
-        n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::double_t>();
       }
     } else if (n.is_type<language::B_set>()) {
       n.m_data_type =
diff --git a/src/language/node_processor/AffectationProcessor.hpp b/src/language/node_processor/AffectationProcessor.hpp
index 59ac02e78d2ac0c77a9138c8e5cd6cb341df35a5..a2c106857f0016953e50a951d655404f514e7712 100644
--- a/src/language/node_processor/AffectationProcessor.hpp
+++ b/src/language/node_processor/AffectationProcessor.hpp
@@ -168,7 +168,93 @@ class AffectationExecutor final : public IAffectationExecutor
 };
 
 template <typename OperatorT, typename ArrayT, typename ValueT, typename DataT>
-class ComponentAffectationExecutor final : public IAffectationExecutor
+class MatrixComponentAffectationExecutor final : public IAffectationExecutor
+{
+ private:
+  ArrayT& m_lhs_array;
+  ASTNode& m_index0_expression;
+  ASTNode& m_index1_expression;
+
+  static inline const bool m_is_defined{[] {
+    if constexpr (not std::is_same_v<typename ArrayT::data_type, ValueT>) {
+      return false;
+    } else if constexpr (std::is_same_v<std::decay_t<ValueT>, bool>) {
+      if constexpr (not std::is_same_v<OperatorT, language::eq_op>) {
+        return false;
+      }
+    }
+    return true;
+  }()};
+
+ public:
+  MatrixComponentAffectationExecutor(ASTNode& node,
+                                     ArrayT& lhs_array,
+                                     ASTNode& index0_expression,
+                                     ASTNode& index1_expression)
+    : m_lhs_array{lhs_array}, m_index0_expression{index0_expression}, m_index1_expression{index1_expression}
+  {
+    // LCOV_EXCL_START
+    if constexpr (not m_is_defined) {
+      throw ParseError("unexpected error: invalid operands to affectation expression", std::vector{node.begin()});
+    }
+    // LCOV_EXCL_STOP
+  }
+
+  PUGS_INLINE void
+  affect(ExecutionPolicy& exec_policy, DataVariant&& rhs)
+  {
+    if constexpr (m_is_defined) {
+      auto get_index_value = [&](DataVariant&& value_variant) -> int64_t {
+        int64_t index_value = 0;
+        std::visit(
+          [&](auto&& value) {
+            using IndexValueT = std::decay_t<decltype(value)>;
+            if constexpr (std::is_integral_v<IndexValueT>) {
+              index_value = value;
+            } else {
+              // LCOV_EXCL_START
+              throw UnexpectedError("invalid index type");
+              // LCOV_EXCL_STOP
+            }
+          },
+          value_variant);
+        return index_value;
+      };
+
+      const int64_t index0_value = get_index_value(m_index0_expression.execute(exec_policy));
+      const int64_t index1_value = get_index_value(m_index1_expression.execute(exec_policy));
+
+      if constexpr (std::is_same_v<ValueT, std::string>) {
+        if constexpr (std::is_same_v<OperatorT, language::eq_op>) {
+          if constexpr (std::is_same_v<std::string, DataT>) {
+            m_lhs_array(index0_value, index1_value) = std::get<DataT>(rhs);
+          } else {
+            m_lhs_array(index0_value, index1_value) = std::to_string(std::get<DataT>(rhs));
+          }
+        } else {
+          if constexpr (std::is_same_v<std::string, DataT>) {
+            m_lhs_array(index0_value, index1_value) += std::get<std::string>(rhs);
+          } else {
+            m_lhs_array(index0_value, index1_value) += std::to_string(std::get<DataT>(rhs));
+          }
+        }
+      } else {
+        if constexpr (std::is_same_v<OperatorT, language::eq_op>) {
+          if constexpr (std::is_same_v<ValueT, DataT>) {
+            m_lhs_array(index0_value, index1_value) = std::get<DataT>(rhs);
+          } else {
+            m_lhs_array(index0_value, index1_value) = static_cast<ValueT>(std::get<DataT>(rhs));
+          }
+        } else {
+          AffOp<OperatorT>().eval(m_lhs_array(index0_value, index1_value), std::get<DataT>(rhs));
+        }
+      }
+    }
+  }
+};
+
+template <typename OperatorT, typename ArrayT, typename ValueT, typename DataT>
+class VectorComponentAffectationExecutor final : public IAffectationExecutor
 {
  private:
   ArrayT& m_lhs_array;
@@ -186,7 +272,7 @@ class ComponentAffectationExecutor final : public IAffectationExecutor
   }()};
 
  public:
-  ComponentAffectationExecutor(ASTNode& node, ArrayT& lhs_array, ASTNode& index_expression)
+  VectorComponentAffectationExecutor(ASTNode& node, ArrayT& lhs_array, ASTNode& index_expression)
     : m_lhs_array{lhs_array}, m_index_expression{index_expression}
   {
     // LCOV_EXCL_START
@@ -289,53 +375,97 @@ class AffectationProcessor final : public INodeProcessor
       Assert(found);
       DataVariant& value = i_symbol->attributes().value();
 
-      // LCOV_EXCL_START
-      if (array_expression.m_data_type != ASTNodeDataType::vector_t) {
-        throw ParseError("unexpected error: invalid lhs (expecting R^d)",
-                         std::vector{array_subscript_expression.begin()});
-      }
-      // LCOV_EXCL_STOP
+      if (array_expression.m_data_type == ASTNodeDataType::vector_t) {
+        Assert(array_subscript_expression.children.size() == 2);
 
-      auto& index_expression = *array_subscript_expression.children[1];
+        auto& index_expression = *array_subscript_expression.children[1];
 
-      switch (array_expression.m_data_type.dimension()) {
-      case 1: {
-        using ArrayTypeT = TinyVector<1>;
-        if (not std::holds_alternative<ArrayTypeT>(value)) {
-          value = ArrayTypeT{};
+        switch (array_expression.m_data_type.dimension()) {
+        case 1: {
+          using ArrayTypeT = TinyVector<1>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor =
+            std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
+          break;
         }
-        using AffectationExecutorT = ComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
-        m_affectation_executor =
-          std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
-        break;
-      }
-      case 2: {
-        using ArrayTypeT = TinyVector<2>;
-        if (not std::holds_alternative<ArrayTypeT>(value)) {
-          value = ArrayTypeT{};
+        case 2: {
+          using ArrayTypeT = TinyVector<2>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor =
+            std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
+          break;
         }
-        using AffectationExecutorT = ComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
-        m_affectation_executor =
-          std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
-        break;
-      }
-      case 3: {
-        using ArrayTypeT = TinyVector<3>;
-        if (not std::holds_alternative<ArrayTypeT>(value)) {
-          value = ArrayTypeT{};
+        case 3: {
+          using ArrayTypeT = TinyVector<3>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor =
+            std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
+          break;
         }
-        using AffectationExecutorT = ComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
-        m_affectation_executor =
-          std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value), index_expression);
-        break;
-      }
-        // LCOV_EXCL_START
-      default: {
-        throw ParseError("unexpected error: invalid vector dimension", std::vector{array_subscript_expression.begin()});
-      }
-        // LCOV_EXCL_STOP
+          // LCOV_EXCL_START
+        default: {
+          throw ParseError("unexpected error: invalid vector dimension",
+                           std::vector{array_subscript_expression.begin()});
+        }
+          // LCOV_EXCL_STOP
+        }
+      } else if (array_expression.m_data_type == ASTNodeDataType::matrix_t) {
+        Assert(array_subscript_expression.children.size() == 3);
+        Assert(array_expression.m_data_type.nbRows() == array_expression.m_data_type.nbColumns());
+
+        auto& index0_expression = *array_subscript_expression.children[1];
+        auto& index1_expression = *array_subscript_expression.children[2];
+
+        switch (array_expression.m_data_type.nbRows()) {
+        case 1: {
+          using ArrayTypeT = TinyMatrix<1>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor     = std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value),
+                                                                          index0_expression, index1_expression);
+          break;
+        }
+        case 2: {
+          using ArrayTypeT = TinyMatrix<2>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor     = std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value),
+                                                                          index0_expression, index1_expression);
+          break;
+        }
+        case 3: {
+          using ArrayTypeT = TinyMatrix<3>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor     = std::make_unique<AffectationExecutorT>(node, std::get<ArrayTypeT>(value),
+                                                                          index0_expression, index1_expression);
+          break;
+        }
+          // LCOV_EXCL_START
+        default: {
+          throw ParseError("unexpected error: invalid vector dimension",
+                           std::vector{array_subscript_expression.begin()});
+        }
+          // LCOV_EXCL_STOP
+        }
+      } else {
+        throw UnexpectedError("");
       }
-
     } else {
       // LCOV_EXCL_START
       throw ParseError("unexpected error: invalid lhs", std::vector{node.children[0]->begin()});
@@ -738,51 +868,96 @@ class ListAffectationProcessor final : public INodeProcessor
       Assert(found);
       DataVariant& value = i_symbol->attributes().value();
 
-      if (array_expression.m_data_type != ASTNodeDataType::vector_t) {
-        // LCOV_EXCL_START
-        throw ParseError("unexpected error: invalid lhs (expecting R^d)",
-                         std::vector{array_subscript_expression.begin()});
-        // LCOV_EXCL_STOP
-      }
+      if (array_expression.m_data_type == ASTNodeDataType::vector_t) {
+        Assert(array_subscript_expression.children.size() == 2);
 
-      auto& index_expression = *array_subscript_expression.children[1];
+        auto& index_expression = *array_subscript_expression.children[1];
 
-      switch (array_expression.m_data_type.dimension()) {
-      case 1: {
-        using ArrayTypeT = TinyVector<1>;
-        if (not std::holds_alternative<ArrayTypeT>(value)) {
-          value = ArrayTypeT{};
+        switch (array_expression.m_data_type.dimension()) {
+        case 1: {
+          using ArrayTypeT = TinyVector<1>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor_list.emplace_back(
+            std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression));
+          break;
         }
-        using AffectationExecutorT = ComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
-        m_affectation_executor_list.emplace_back(
-          std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression));
-        break;
-      }
-      case 2: {
-        using ArrayTypeT = TinyVector<2>;
-        if (not std::holds_alternative<ArrayTypeT>(value)) {
-          value = ArrayTypeT{};
+        case 2: {
+          using ArrayTypeT = TinyVector<2>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor_list.emplace_back(
+            std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression));
+          break;
         }
-        using AffectationExecutorT = ComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
-        m_affectation_executor_list.emplace_back(
-          std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression));
-        break;
-      }
-      case 3: {
-        using ArrayTypeT = TinyVector<3>;
-        if (not std::holds_alternative<ArrayTypeT>(value)) {
-          value = ArrayTypeT{};
+        case 3: {
+          using ArrayTypeT = TinyVector<3>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = VectorComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor_list.emplace_back(
+            std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression));
+          break;
+        }
+          // LCOV_EXCL_START
+        default: {
+          throw ParseError("unexpected error: invalid vector dimension",
+                           std::vector{array_subscript_expression.begin()});
+        }
+          // LCOV_EXCL_STOP
+        }
+      } else if (array_expression.m_data_type == ASTNodeDataType::matrix_t) {
+        Assert(array_subscript_expression.children.size() == 3);
+
+        auto& index0_expression = *array_subscript_expression.children[1];
+        auto& index1_expression = *array_subscript_expression.children[2];
+
+        switch (array_expression.m_data_type.dimension()) {
+        case 1: {
+          using ArrayTypeT = TinyMatrix<1>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor_list.emplace_back(
+            std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index0_expression,
+                                                   index1_expression));
+          break;
+        }
+        case 2: {
+          using ArrayTypeT = TinyMatrix<2>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor_list.emplace_back(
+            std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index0_expression,
+                                                   index1_expression));
+          break;
+        }
+        case 3: {
+          using ArrayTypeT = TinyMatrix<3>;
+          if (not std::holds_alternative<ArrayTypeT>(value)) {
+            value = ArrayTypeT{};
+          }
+          using AffectationExecutorT = MatrixComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
+          m_affectation_executor_list.emplace_back(
+            std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index0_expression,
+                                                   index1_expression));
+          break;
+        }
+          // LCOV_EXCL_START
+        default: {
+          throw ParseError("unexpected error: invalid vector dimension",
+                           std::vector{array_subscript_expression.begin()});
+        }
+          // LCOV_EXCL_STOP
         }
-        using AffectationExecutorT = ComponentAffectationExecutor<OperatorT, ArrayTypeT, ValueT, DataT>;
-        m_affectation_executor_list.emplace_back(
-          std::make_unique<AffectationExecutorT>(lhs_node, std::get<ArrayTypeT>(value), index_expression));
-        break;
-      }
-        // LCOV_EXCL_START
-      default: {
-        throw ParseError("unexpected error: invalid vector dimension", std::vector{array_subscript_expression.begin()});
-      }
-        // LCOV_EXCL_STOP
       }
     } else {
       // LCOV_EXCL_START
diff --git a/src/language/node_processor/ArraySubscriptProcessor.hpp b/src/language/node_processor/ArraySubscriptProcessor.hpp
index 3a71c01846b9e17ce28818d732716aea764d2496..43eb962b9d724b19b4cfe7ab808f52dd73d1fce2 100644
--- a/src/language/node_processor/ArraySubscriptProcessor.hpp
+++ b/src/language/node_processor/ArraySubscriptProcessor.hpp
@@ -15,9 +15,7 @@ class ArraySubscriptProcessor : public INodeProcessor
   DataVariant
   execute(ExecutionPolicy& exec_policy)
   {
-    auto& index_expression = *m_array_subscript_expression.children[1];
-
-    const int64_t index_value = [&](DataVariant&& value_variant) -> int64_t {
+    auto get_index_value = [&](DataVariant&& value_variant) -> int64_t {
       int64_t index_value = 0;
       std::visit(
         [&](auto&& value) {
@@ -26,20 +24,39 @@ class ArraySubscriptProcessor : public INodeProcessor
             index_value = value;
           } else {
             // LCOV_EXCL_START
-            throw ParseError("unexpected error: invalid index type", std::vector{index_expression.begin()});
+            throw UnexpectedError("invalid index type");
             // LCOV_EXCL_STOP
           }
         },
         value_variant);
       return index_value;
-    }(index_expression.execute(exec_policy));
+    };
+
+    if constexpr (is_tiny_vector_v<ArrayTypeT>) {
+      auto& index_expression = *m_array_subscript_expression.children[1];
+
+      const int64_t index_value = get_index_value(index_expression.execute(exec_policy));
+
+      auto& array_expression = *m_array_subscript_expression.children[0];
+
+      auto&& array_value = array_expression.execute(exec_policy);
+      ArrayTypeT& array  = std::get<ArrayTypeT>(array_value);
+
+      return array[index_value];
+    } else if constexpr (is_tiny_matrix_v<ArrayTypeT>) {
+      auto& index0_expression = *m_array_subscript_expression.children[1];
+      auto& index1_expression = *m_array_subscript_expression.children[2];
+
+      const int64_t index0_value = get_index_value(index0_expression.execute(exec_policy));
+      const int64_t index1_value = get_index_value(index1_expression.execute(exec_policy));
 
-    auto& array_expression = *m_array_subscript_expression.children[0];
+      auto& array_expression = *m_array_subscript_expression.children[0];
 
-    auto&& array_value = array_expression.execute(exec_policy);
-    ArrayTypeT& array  = std::get<ArrayTypeT>(array_value);
+      auto&& array_value = array_expression.execute(exec_policy);
+      ArrayTypeT& array  = std::get<ArrayTypeT>(array_value);
 
-    return array[index_value];
+      return array(index0_value, index1_value);
+    }
   }
 
   ArraySubscriptProcessor(ASTNode& array_subscript_expression)
diff --git a/tests/test_ASTNodeDataTypeBuilder.cpp b/tests/test_ASTNodeDataTypeBuilder.cpp
index 6c37b1afe9a7c411e2c1643ac1bd1e3cb1a72aeb..f12e31ab4f9378777b2c4da18297aeb51e932ff8 100644
--- a/tests/test_ASTNodeDataTypeBuilder.cpp
+++ b/tests/test_ASTNodeDataTypeBuilder.cpp
@@ -305,7 +305,7 @@ let x : R; x[2];
         auto ast = ASTBuilder::build(input);
         ASTSymbolTableBuilder{*ast};
 
-        REQUIRE_THROWS_WITH(ASTNodeDataTypeBuilder{*ast}, "invalid types 'R[Z]' for array subscript");
+        REQUIRE_THROWS_WITH(ASTNodeDataTypeBuilder{*ast}, "invalid subscript expression: R cannot be indexed");
       }
 
       SECTION("too many variables")