diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp
index 409c4427d..7a9d8867c 100644
--- a/src/solver/EigenLinearEquationSolver.cpp
+++ b/src/solver/EigenLinearEquationSolver.cpp
@@ -130,6 +130,12 @@ namespace storm {
             // Typedef the map-type so we don't have to spell it out.
             typedef decltype(Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
             
+            bool multiplyResultProvided = multiplyResult != nullptr;
+            if (!multiplyResult) {
+                multiplyResult = new std::vector<ValueType>(eigenA->cols());
+            }
+            auto eigenMultiplyResult = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(multiplyResult->data(), multiplyResult->size());
+            
             // Map the input vectors x and b to Eigen's format.
             std::unique_ptr<MapType> eigenB;
             if (b != nullptr) {
@@ -138,11 +144,23 @@ namespace storm {
             auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             
             // Perform n matrix-vector multiplications.
+            auto currentX = &eigenX;
+            auto nextX = &eigenMultiplyResult;
             for (uint64_t iteration = 0; iteration < n; ++iteration) {
-                eigenX = *this->eigenA * eigenX;
+                *nextX = *eigenA * *currentX;
                 if (eigenB != nullptr) {
-                    eigenX += *eigenB;
+                    *nextX += *eigenB;
                 }
+                std::swap(nextX, currentX);
+            }
+            
+            // If the last result we obtained is not the one in the input vector x, we swap the result there.
+            if (currentX != &eigenX) {
+                std::swap(*nextX, *currentX);
+            }
+            
+            if (!multiplyResultProvided) {
+                delete multiplyResult;
             }
         }
         
@@ -174,19 +192,37 @@ namespace storm {
             // Typedef the map-type so we don't have to spell it out.
             typedef decltype(Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
 
+            bool multiplyResultProvided = multiplyResult != nullptr;
+            if (!multiplyResult) {
+                multiplyResult = new std::vector<storm::RationalNumber>(eigenA->cols());
+            }
+            auto eigenMultiplyResult = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(multiplyResult->data(), multiplyResult->size());
+            
             // Map the input vectors x and b to Eigen's format.
             std::unique_ptr<MapType> eigenB;
             if (b != nullptr) {
                 eigenB = std::make_unique<MapType>(Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b->data(), b->size()));
             }
-            Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
+            auto eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             
             // Perform n matrix-vector multiplications.
+            auto currentX = &eigenX;
+            auto nextX = &eigenMultiplyResult;
             for (uint64_t iteration = 0; iteration < n; ++iteration) {
-                eigenX = *eigenA * eigenX;
+                *nextX = *eigenA * *currentX;
                 if (eigenB != nullptr) {
-                    eigenX += *eigenB;
+                    *nextX += *eigenB;
                 }
+                std::swap(currentX, nextX);
+            }
+            
+            // If the last result we obtained is not the one in the input vector x, we swap the result there.
+            if (currentX != &eigenX) {
+                std::swap(*nextX, *currentX);
+            }
+            
+            if (!multiplyResultProvided) {
+                delete multiplyResult;
             }
         }
         
@@ -208,19 +244,37 @@ namespace storm {
             // Typedef the map-type so we don't have to spell it out.
             typedef decltype(Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
 
+            bool multiplyResultProvided = multiplyResult != nullptr;
+            if (!multiplyResult) {
+                multiplyResult = new std::vector<storm::RationalFunction>(eigenA->cols());
+            }
+            auto eigenMultiplyResult = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(multiplyResult->data(), multiplyResult->size());
+            
             // Map the input vectors x and b to Eigen's format.
             std::unique_ptr<MapType> eigenB;
             if (b != nullptr) {
                 eigenB = std::make_unique<MapType>(Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b->data(), b->size()));
             }
-            Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
+            auto eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             
             // Perform n matrix-vector multiplications.
+            auto currentX = &eigenX;
+            auto nextX = &eigenMultiplyResult;
             for (uint64_t iteration = 0; iteration < n; ++iteration) {
-                eigenX = *eigenA * eigenX;
+                *nextX = *eigenA * *currentX;
                 if (eigenB != nullptr) {
-                    eigenX += *eigenB;
+                    *nextX += *eigenB;
                 }
+                std::swap(nextX, currentX);
+            }
+            
+            // If the last result we obtained is not the one in the input vector x, we swap the result there.
+            if (currentX != &eigenX) {
+                std::swap(*nextX, *currentX);
+            }
+            
+            if (!multiplyResultProvided) {
+                delete multiplyResult;
             }
         }
         
diff --git a/test/functional/builder/parametric_die.pm b/test/functional/builder/parametric_die.pm
new file mode 100644
index 000000000..bd4e8b103
--- /dev/null
+++ b/test/functional/builder/parametric_die.pm
@@ -0,0 +1,34 @@
+// Knuth's model of a fair die using only fair coins
+dtmc
+
+const double p;
+
+module die
+
+	// local state
+	s : [0..7] init 0;
+	// value of the dice
+	d : [0..6] init 0;
+	
+	[] s=0 -> p : (s'=1) + (1-p) : (s'=2);
+	[] s=1 -> p : (s'=3) + (1-p) : (s'=4);
+	[] s=2 -> p : (s'=5) + (1-p) : (s'=6);
+	[] s=3 -> p : (s'=1) + (1-p) : (s'=7) & (d'=1);
+	[] s=4 -> p : (s'=7) & (d'=2) + (1-p) : (s'=7) & (d'=3);
+	[] s=5 -> p : (s'=7) & (d'=4) + (1-p) : (s'=7) & (d'=5);
+	[] s=6 -> p : (s'=2) + (1-p) : (s'=7) & (d'=6);
+	[] s=7 -> 1: (s'=7);
+	
+endmodule
+
+rewards "coin_flips"
+	[] s<7 : 1;
+endrewards
+
+label "one" = s=7&d=1;
+label "two" = s=7&d=2;
+label "three" = s=7&d=3;
+label "four" = s=7&d=4;
+label "five" = s=7&d=5;
+label "six" = s=7&d=6;
+label "done" = s=7;
diff --git a/test/functional/solver/EigenLinearEquationSolverTest.cpp b/test/functional/solver/EigenLinearEquationSolverTest.cpp
index 2bf575c2a..5bc1b4a42 100644
--- a/test/functional/solver/EigenLinearEquationSolverTest.cpp
+++ b/test/functional/solver/EigenLinearEquationSolverTest.cpp
@@ -207,7 +207,7 @@ TEST(DISABLED_EigenLinearEquationSolver, BiCGSTAB_Diagonal) {
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
 }
 
-TEST(EigenLinearEquationSolver, MatrixVectorMultplication) {
+TEST(EigenLinearEquationSolver, MatrixVectorMultiplication) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));
diff --git a/test/functional/solver/EliminationLinearEquationSolverTest.cpp b/test/functional/solver/EliminationLinearEquationSolverTest.cpp
index 39007a5d7..7a964a1d6 100644
--- a/test/functional/solver/EliminationLinearEquationSolverTest.cpp
+++ b/test/functional/solver/EliminationLinearEquationSolverTest.cpp
@@ -34,7 +34,7 @@ TEST(EliminationLinearEquationSolver, Solve) {
     ASSERT_LT(std::abs(x[2] - (-1)), 1e-15);
 }
 
-TEST(EliminationLinearEquationSolver, MatrixVectorMultplication) {
+TEST(EliminationLinearEquationSolver, MatrixVectorMultiplication) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));
diff --git a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
index 181d3429a..44282f522 100644
--- a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
+++ b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
@@ -227,7 +227,7 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) {
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
 }
 
-TEST(GmmxxLinearEquationSolver, MatrixVectorMultplication) {
+TEST(GmmxxLinearEquationSolver, MatrixVectorMultiplication) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));
diff --git a/test/functional/solver/NativeLinearEquationSolverTest.cpp b/test/functional/solver/NativeLinearEquationSolverTest.cpp
index a1313a797..cfe61741c 100644
--- a/test/functional/solver/NativeLinearEquationSolverTest.cpp
+++ b/test/functional/solver/NativeLinearEquationSolverTest.cpp
@@ -34,7 +34,7 @@ TEST(NativeLinearEquationSolver, SolveWithStandardOptions) {
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
 }
 
-TEST(NativeLinearEquationSolver, MatrixVectorMultplication) {
+TEST(NativeLinearEquationSolver, MatrixVectorMultiplication) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));