diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
index 0a6dabcd8..3b6e4fb70 100644
--- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
@@ -205,7 +205,12 @@ namespace storm {
             std::vector<ValueType> SparseCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) {
                 return SparseDtmcPrctlHelper<ValueType>::computeUntilProbabilities(env, std::move(goal), computeProbabilityMatrix(rateMatrix, exitRateVector), backwardTransitions, phiStates, psiStates, qualitative);
             }
-            
+
+            template <typename ValueType>
+            std::vector<ValueType> SparseCtmcCslHelper::computeAllUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                return SparseDtmcPrctlHelper<ValueType>::computeAllUntilProbabilities(env, std::move(goal), computeProbabilityMatrix(rateMatrix, exitRateVector), backwardTransitions, initialStates, phiStates, psiStates);
+            }
+
             template <typename ValueType>
             std::vector<ValueType> SparseCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& nextStates) {
                 return SparseDtmcPrctlHelper<ValueType>::computeNextProbabilities(env, computeProbabilityMatrix(rateMatrix, exitRateVector), nextStates);
@@ -808,6 +813,8 @@ namespace storm {
             
             template std::vector<double> SparseCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, storm::storage::SparseMatrix<double> const& backwardTransitions, std::vector<double> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative);
 
+            template std::vector<double> SparseCtmcCslHelper::computeAllUntilProbabilities(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, storm::storage::SparseMatrix<double> const& backwardTransitions, std::vector<double> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates);
+
             template std::vector<double> SparseCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& exitRateVector, storm::storage::BitVector const& nextStates);
             
             template std::vector<double> SparseCtmcCslHelper::computeInstantaneousRewards(Environment const& env, storm::solver::SolveGoal<double>&& goal, storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, double timeBound);
@@ -835,6 +842,9 @@ namespace storm {
             template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalNumber> const& backwardTransitions, std::vector<storm::RationalNumber> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative);
             template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalFunction> const& backwardTransitions, std::vector<storm::RationalFunction> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative);
 
+            template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeAllUntilProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalNumber>&& goal, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalNumber> const& backwardTransitions, std::vector<storm::RationalNumber> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates);
+            template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeAllUntilProbabilities(Environment const& env, storm::solver::SolveGoal<storm::RationalFunction>&& goal, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalFunction> const& backwardTransitions, std::vector<storm::RationalFunction> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates);
+
             template std::vector<storm::RationalNumber> SparseCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::storage::SparseMatrix<storm::RationalNumber> const& rateMatrix, std::vector<storm::RationalNumber> const& exitRateVector, storm::storage::BitVector const& nextStates);
             template std::vector<storm::RationalFunction> SparseCtmcCslHelper::computeNextProbabilities(Environment const& env, storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, std::vector<storm::RationalFunction> const& exitRateVector, storm::storage::BitVector const& nextStates);
 
diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h
index 4ecd66abc..233e695ff 100644
--- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h
+++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.h
@@ -27,6 +27,9 @@ namespace storm {
                 template <typename ValueType>
                 static std::vector<ValueType> computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative);
 
+                template <typename ValueType>
+                static std::vector<ValueType> computeAllUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates);
+
                 template <typename ValueType>
                 static std::vector<ValueType> computeNextProbabilities(Environment const& env, storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& nextStates);
                 
diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
index 03d3fb1a0..dd8bbce8f 100644
--- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
@@ -330,7 +330,87 @@ namespace storm {
                 }
                 return result;
             }
-            
+
+            template<typename ValueType, typename RewardModelType>
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeAllUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+
+                uint_fast64_t numberOfStates = transitionMatrix.getRowCount();
+                std::vector<ValueType> result(numberOfStates, storm::utility::zero<ValueType>());
+
+                // We need to identify the maybe states (states which have a probability for satisfying the until formula
+                // that is strictly between 0 and 1) and the states that satisfy the formula with probability 1.
+                //storm::storage::BitVector maybeStates, statesWithProbability1;
+                storm::storage::BitVector relevantStates(numberOfStates, true);
+
+                // Get all states that have probability 0 and 1 of satisfying the until-formula.
+                //std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(transitionMatrix, phiStates, initialStates);
+                //storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
+                //statesWithProbability1 = std::move(statesWithProbability01.second);
+                //maybeStates = ~(statesWithProbability0 | statesWithProbability1);
+
+                //STORM_LOG_INFO("Preprocessing: " << statesWithProbability1.getNumberOfSetBits() << " states with probability 1, " << statesWithProbability0.getNumberOfSetBits() << " with probability 0 (" << maybeStates.getNumberOfSetBits() << " states remaining).");
+
+                // Set values of resulting vector that are known exactly.
+                //storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability0, storm::utility::zero<ValueType>());
+                //storm::utility::vector::setVectorValues<ValueType>(result, statesWithProbability1, storm::utility::one<ValueType>());
+
+                // Compute exact probabilities for some states.
+                if (!relevantStates.empty()) {
+                    // In this case we have to compute the probabilities.
+
+                    // Check whether we need to convert the input to equation system format.
+                    storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
+                    bool convertToEquationSystem = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem;
+
+                    // We can eliminate the rows and columns from the original transition probability matrix.
+                    storm::storage::SparseMatrix<ValueType> submatrix(transitionMatrix);
+                    submatrix.makeRowsAbsorbing(psiStates);
+                    submatrix.deleteDiagonalEntries(psiStates);
+                    storm::storage::BitVector failState(numberOfStates, false);
+                    failState.set(0, true);
+                    submatrix.deleteDiagonalEntries(failState);
+                    submatrix = submatrix.transpose();
+                    submatrix = submatrix.getSubmatrix(true, relevantStates, relevantStates, convertToEquationSystem);
+
+                    if (convertToEquationSystem) {
+                        // Converting the matrix from the fixpoint notation to the form needed for the equation
+                        // system. That is, we go from x = A*x + b to (I-A)x = b.
+                        submatrix.convertToEquationSystem();
+                    }
+
+                    // Initialize the x vector with 0.5 for each element.
+                    // This is the initial guess for the iterative solvers. It should be safe as for all
+                    // 'maybe' states we know that the probability is strictly larger than 0.
+                    std::vector<ValueType> x = std::vector<ValueType>(relevantStates.getNumberOfSetBits(), storm::utility::convertNumber<ValueType>(0.5));
+                    //std::vector<ValueType> x = std::vector<ValueType>(relevantStates.getNumberOfSetBits(), storm::utility::zero<ValueType>());
+
+                    std::vector<ValueType> b(relevantStates.getNumberOfSetBits(), storm::utility::zero<ValueType>());
+                    // Set initial states
+                    size_t i = 0;
+                    ValueType initDist = storm::utility::one<ValueType>() / storm::utility::convertNumber<ValueType>(initialStates.getNumberOfSetBits());
+                    for (auto const& state : relevantStates) {
+                        if (initialStates.get(state)) {
+                            b[i] = initDist;
+                        }
+                        ++i;
+                    }
+
+                    // Prepare the right-hand side of the equation system. For entry i this corresponds to
+                    // the accumulated probability of going from state i to some 'yes' state.
+                    //std::vector<ValueType> b = transitionMatrix.getConstrainedRowSumVector(relevantStates, statesWithProbability1);
+
+                    // Now solve the created system of linear equations.
+                    goal.restrictRelevantValues(relevantStates);
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = storm::solver::configureLinearEquationSolver(env, std::move(goal), linearEquationSolverFactory, std::move(submatrix));
+                    solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
+                    solver->solveEquations(env, x, b);
+
+                    // Set values of resulting vector according to result.
+                    storm::utility::vector::setVectorValues<ValueType>(result, relevantStates, x);
+                }
+                return result;
+            }
+
             template<typename ValueType, typename RewardModelType>
             std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeGloballyProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative) {
                 goal.oneMinus();
diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
index a4e2882bc..1492b29df 100644
--- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
+++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
@@ -36,6 +36,8 @@ namespace storm {
                 
                 static std::vector<ValueType> computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, ModelCheckerHint const& hint = ModelCheckerHint());
 
+                static std::vector<ValueType> computeAllUntilProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates);
+
                 static std::vector<ValueType> computeGloballyProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative);
                 
                 static std::vector<ValueType> computeCumulativeRewards(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound);
diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp
index 729b85b46..49cef9e4e 100644
--- a/src/storm/storage/SparseMatrix.cpp
+++ b/src/storm/storage/SparseMatrix.cpp
@@ -1301,6 +1301,19 @@ namespace storm {
                 }
             }
         }
+
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::deleteDiagonalEntries(storm::storage::BitVector const& states) {
+            // Iterate over all rows and negate all the elements that are not on the diagonal.
+            for (index_type group = 0; group < this->getRowGroupCount(); ++group) {
+                for (auto& entry : this->getRowGroup(group)) {
+                    if (entry.getColumn() == group && states[group]) {
+                        --this->nonzeroEntryCount;
+                        entry.setValue(storm::utility::zero<ValueType>());
+                    }
+                }
+            }
+        }
         
         template<typename ValueType>
         typename std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> SparseMatrix<ValueType>::getJacobiDecomposition() const {
diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h
index 4013a7cbf..e61d34624 100644
--- a/src/storm/storage/SparseMatrix.h
+++ b/src/storm/storage/SparseMatrix.h
@@ -799,6 +799,11 @@ namespace storm {
              * Sets all diagonal elements to zero.
              */
             void deleteDiagonalEntries();
+
+            /*!
+             * Sets all diagonal elements to zero.
+             */
+            void deleteDiagonalEntries(storm::storage::BitVector const& states);
             
             /*!
              * Calculates the Jacobi decomposition of this sparse matrix. For this operation, the matrix must be square.