diff --git a/3rdparty/simlib/include/mata/simlib/explicit_lts.hh b/3rdparty/simlib/include/mata/simlib/explicit_lts.hh
index 94cbe2c6..c5aaccf1 100644
--- a/3rdparty/simlib/include/mata/simlib/explicit_lts.hh
+++ b/3rdparty/simlib/include/mata/simlib/explicit_lts.hh
@@ -21,6 +21,9 @@ namespace Simlib { class ExplicitLTS; }
 
 class Simlib::ExplicitLTS {
 
+    /// Map mapping label to a unique number in interval [0, <number of labels in automaton>].
+    /// The unique number is used as an index for accessing data of explicit LTS corresponding to the given label.
+    std::unordered_map<size_t, size_t> symbol_map{};
 	size_t states_{};
 	size_t transitions_;
 	std::vector<
diff --git a/3rdparty/simlib/src/explicit_lts_sim.cc b/3rdparty/simlib/src/explicit_lts_sim.cc
index 6a081f3c..cd45603c 100644
--- a/3rdparty/simlib/src/explicit_lts_sim.cc
+++ b/3rdparty/simlib/src/explicit_lts_sim.cc
@@ -903,11 +903,9 @@ BinaryRelation Simlib::ExplicitLTS::compute_simulation()
 }
 
 
-void Simlib::ExplicitLTS::add_transition(
-	size_t   q,
-	size_t   a,
-	size_t   r)
-{
+void Simlib::ExplicitLTS::add_transition(size_t q, size_t a, size_t r) {
+	if (!symbol_map.contains(a)) { symbol_map[a] = symbol_map.size(); };
+	a = symbol_map[a];
 	if (a >= this->data_.size())
 	{
 		this->data_.resize(a + 1);
diff --git a/src/nfa/operations.cc b/src/nfa/operations.cc
index a0aabbd8..2d25e411 100644
--- a/src/nfa/operations.cc
+++ b/src/nfa/operations.cc
@@ -38,15 +38,12 @@ using StateBoolArray = std::vector<bool>; ///< Bool array for states in the auto
 
 namespace {
     Simlib::Util::BinaryRelation compute_fw_direct_simulation(const Nfa& aut) {
-        Symbol maxSymbol = 0;
+        Symbol maxSymbol{ aut.get_max_symbol() };
         const size_t state_num = aut.size();
         Simlib::ExplicitLTS LTSforSimulation(state_num);
 
         for (const Transition& transition : aut.delta.transitions()) {
             LTSforSimulation.add_transition(transition.source, transition.symbol, transition.target);
-            if (transition.symbol > maxSymbol) {
-                maxSymbol = transition.symbol;
-            }
         }
 
         // final states cannot be simulated by nonfinal -> we add new selfloops over final states with new symbol in LTS
diff --git a/tests-integration/src/reduce_epsilon.cc b/tests-integration/src/reduce_epsilon.cc
new file mode 100644
index 00000000..5d82eb0f
--- /dev/null
+++ b/tests-integration/src/reduce_epsilon.cc
@@ -0,0 +1,39 @@
+#include "utils/utils.hh"
+
+#include "mata/nfa/nfa.hh"
+
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <chrono>
+#include <string>
+#include <cstring>
+
+using namespace mata::nfa;
+
+int main() {
+    Nfa b{10};
+    b.initial.insert(0);
+    b.final.insert({2, 4, 8, 7});
+    b.delta.add(0, 'b', 1);
+    b.delta.add(0, 'a', 2);
+    b.delta.add(2, 'a', 4);
+    b.delta.add(2, EPSILON, 3);
+    b.delta.add(3, 'b', 4);
+    b.delta.add(0, 'c', 5);
+    b.delta.add(5, 'a', 8);
+    b.delta.add(5, EPSILON, 6);
+    b.delta.add(6, 'a', 9);
+    b.delta.add(6, 'b', 7);
+
+    Nfa c;
+    TIME_BLOCK(reduce,
+        for (size_t i{ 0 }; i < 1'000; ++i) {
+            c = reduce(b);
+        }
+    );
+    if (!are_equivalent(b, c)) { // Check correctness.
+        return EXIT_FAILURE;
+    }
+}