diff --git a/src/modelchecker/GmmxxMdpPrctlModelChecker.h b/src/modelchecker/GmmxxMdpPrctlModelChecker.h
index 90526a24c..fce08e1c7 100644
--- a/src/modelchecker/GmmxxMdpPrctlModelChecker.h
+++ b/src/modelchecker/GmmxxMdpPrctlModelChecker.h
@@ -59,16 +59,12 @@ private:
 	 * @param A The matrix that is to be multiplied against the vector.
 	 * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter,
 	 * i.e. after the method returns, this vector will contain the computed values.
+     * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix.
 	 * @param b If not null, this vector is being added to the result after each matrix-vector multiplication.
 	 * @param n Specifies the number of iterations the matrix-vector multiplication is performed.
 	 * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector.
 	 */
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const {
-		// Get the starting row indices for the non-deterministic choices to reduce the resulting
-		// vector properly.
-		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = *this->getModel().getNondeterministicChoiceIndices();
-
-
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const {
 		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
 		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(A);
 
diff --git a/src/modelchecker/SparseDtmcPrctlModelChecker.h b/src/modelchecker/SparseDtmcPrctlModelChecker.h
index feb77a88d..d574613e4 100644
--- a/src/modelchecker/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/SparseDtmcPrctlModelChecker.h
@@ -90,17 +90,20 @@ public:
 
         // If we identify the states that have probability 0 of reaching the target states, we can exclude them in the
         // further analysis.
-        storm::storage::BitVector maybeStates = storm::utility::graph::performProbGreater0(this->getModel(), *leftStates, *rightStates);
+        storm::storage::BitVector maybeStates = storm::utility::graph::performProbGreater0(this->getModel(), *leftStates, *rightStates, true, formula.getBound());
         
         // Now we can eliminate the rows and columns from the original transition probability matrix that have probability 0.
         storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates);
         
+        // Compute the new set of target states in the reduced system.
+        storm::storage::BitVector rightStatesInReducedSystem = maybeStates % *rightStates;
+        
 		// Make all rows absorbing that satisfy the second sub-formula.
-		submatrix.makeRowsAbsorbing(maybeStates % *rightStates);
+		submatrix.makeRowsAbsorbing(rightStatesInReducedSystem);
 
 		// Create the vector with which to multiply.
         std::vector<Type> subresult(maybeStates.getNumberOfSetBits());
-		storm::utility::vector::setVectorValues(subresult, *rightStates, storm::utility::constGetOne<Type>());
+		storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::constGetOne<Type>());
 
 		// Perform the matrix vector multiplication as often as required by the formula bound.
 		this->performMatrixVectorMultiplication(submatrix, subresult, nullptr, formula.getBound());
diff --git a/src/modelchecker/SparseMdpPrctlModelChecker.h b/src/modelchecker/SparseMdpPrctlModelChecker.h
index 6c1694276..61e3494c6 100644
--- a/src/modelchecker/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/SparseMdpPrctlModelChecker.h
@@ -93,18 +93,37 @@ public:
 		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
 		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
 		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
+        
+        // Determine the states that have 0 probability of reaching the target states.
+        storm::storage::BitVector maybeStates;
+        if (this->minimumOperatorStack.top()) {
+			maybeStates = storm::utility::graph::performProbGreater0A(this->getModel(), this->getModel().getBackwardTransitions(), *leftStates, *rightStates, true, formula.getBound());
+		} else {
+			maybeStates = storm::utility::graph::performProbGreater0E(this->getModel(), this->getModel().getBackwardTransitions(), *leftStates, *rightStates, true, formula.getBound());
+		}
+        
+        // Now we can eliminate the rows and columns from the original transition probability matrix that have probability 0.
+        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices());
 
-		// Copy the matrix before we make any changes.
-		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionMatrix());
+        // Get the "new" nondeterministic choice indices for the submatrix.
+        std::vector<uint_fast64_t> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
 
-		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
-		tmpMatrix.makeRowsAbsorbing(~(*leftStates | *rightStates) | *rightStates, *this->getModel().getNondeterministicChoiceIndices());
-        
-		// Create the vector with which to multiply.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-		storm::utility::vector::setVectorValues(*result, *rightStates, storm::utility::constGetOne<Type>());
+        // Compute the new set of target states in the reduced system.
+        storm::storage::BitVector rightStatesInReducedSystem = maybeStates % *rightStates;
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, nullptr, formula.getBound());
+		// Make all rows absorbing that satisfy the second sub-formula.
+		submatrix.makeRowsAbsorbing(rightStatesInReducedSystem, subNondeterministicChoiceIndices);
+
+        // Create the vector with which to multiply.
+        std::vector<Type> subresult(maybeStates.getNumberOfSetBits());
+		storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::constGetOne<Type>());
+
+		this->performMatrixVectorMultiplication(submatrix, subresult, subNondeterministicChoiceIndices, nullptr, formula.getBound());
+
+		// Create the resulting vector.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+        storm::utility::vector::setVectorValues(*result, maybeStates, subresult);
+		storm::utility::vector::setVectorValues(*result, ~maybeStates, storm::utility::constGetZero<Type>());
 
 		// Delete intermediate results and return result.
 		delete leftStates;
@@ -134,7 +153,7 @@ public:
 		// Delete obsolete sub-result.
 		delete nextStates;
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result);
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, *this->getModel().getNondeterministicChoiceIndices());
 
 		// Return result.
 		return result;
@@ -154,7 +173,7 @@ public:
 	virtual std::vector<Type>* checkBoundedEventually(const storm::property::prctl::BoundedEventually<Type>& formula, bool qualitative) const {
 		// Create equivalent temporary bounded until formula and check it.
 		storm::property::prctl::BoundedUntil<Type> temporaryBoundedUntilFormula(new storm::property::prctl::Ap<Type>("true"), formula.getChild().clone(), formula.getBound());
-		return this->checkBoundedUntil(temporaryBoundedUntilFormula, qualitative);
+        return this->checkBoundedUntil(temporaryBoundedUntilFormula, qualitative);
 	}
 
 	/*!
@@ -286,7 +305,7 @@ public:
 		// Initialize result to state rewards of the model.
 		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, nullptr, formula.getBound());
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, *this->getModel().getNondeterministicChoiceIndices(), nullptr, formula.getBound());
 
 		// Return result.
 		return result;
@@ -329,7 +348,7 @@ public:
 			result = new std::vector<Type>(this->getModel().getNumberOfStates());
 		}
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, &totalRewardVector, formula.getBound());
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, *this->getModel().getNondeterministicChoiceIndices(), &totalRewardVector, formula.getBound());
 
 		// Delete temporary variables and return result.
 		return result;
@@ -446,15 +465,12 @@ private:
 	 * @param A The matrix that is to be multiplied against the vector.
 	 * @param x The initial vector that is to be multiplied against the matrix. This is also the output parameter,
 	 * i.e. after the method returns, this vector will contain the computed values.
+     * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix.
 	 * @param b If not null, this vector is being added to the result after each matrix-vector multiplication.
 	 * @param n Specifies the number of iterations the matrix-vector multiplication is performed.
 	 * @returns The result of the repeated matrix-vector multiplication as the content of the parameter vector.
 	 */
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const {
-		// Get the starting row indices for the non-deterministic choices to reduce the resulting
-		// vector properly.
-		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = *this->getModel().getNondeterministicChoiceIndices();
-
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type>* b = nullptr, uint_fast64_t n = 1) const {
 		// Create vector for result of multiplication, which is reduced to the result vector after
 		// each multiplication.
 		std::vector<Type> multiplyResult(A.getRowCount());
@@ -485,6 +501,7 @@ private:
 	 * @param x The solution vector x. The initial values of x represent a guess of the real values to the solver, but
 	 * may be ignored.
 	 * @param b The right-hand side of the equation system.
+     * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix.
 	 * @returns The solution vector x of the system of linear equations as the content of the parameter x.
 	 */
 	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h
index 7295f29bf..d412b2eb7 100644
--- a/src/storage/BitVector.h
+++ b/src/storage/BitVector.h
@@ -111,8 +111,8 @@ public:
 	BitVector(uint_fast64_t length, bool initTrue = false) : bitCount(length), endIterator(*this, length, length, false), truncateMask((1ll << (bitCount & mod64mask)) - 1ll) {
 		// Check whether the given length is valid.
 		if (length == 0) {
-			LOG4CPLUS_ERROR(logger, "Trying to create bit vector of size 0.");
-			throw storm::exceptions::InvalidArgumentException("Trying to create a bit vector of size 0.");
+			LOG4CPLUS_ERROR(logger, "Cannot create bit vector of size 0.");
+			throw storm::exceptions::InvalidArgumentException() << "Cannot create bit vector of size 0.";
 		}
 
 		// Compute the correct number of buckets needed to store the given number of bits
@@ -272,7 +272,7 @@ public:
 	 */
 	void set(const uint_fast64_t index, const bool value) {
 		uint_fast64_t bucket = index >> 6;
-		if (bucket >= this->bucketCount) throw storm::exceptions::OutOfRangeException();
+		if (bucket >= this->bucketCount) throw storm::exceptions::OutOfRangeException() << "Written index " << index << " out of bounds.";
 		uint_fast64_t mask = static_cast<uint_fast64_t>(1) << (index & mod64mask);
 		if (value) {
 			this->bucketArray[bucket] |= mask;
@@ -290,7 +290,7 @@ public:
 	 */
 	bool get(const uint_fast64_t index) const {
 		uint_fast64_t bucket = index >> 6;
-		if (bucket >= this->bucketCount) throw storm::exceptions::OutOfRangeException();
+		if (bucket >= this->bucketCount) throw storm::exceptions::OutOfRangeException() << "Read index " << index << " out of bounds.";
 		uint_fast64_t mask = static_cast<uint_fast64_t>(1) << (index & mod64mask);
 		return ((this->bucketArray[bucket] & mask) == mask);
 	}
@@ -600,7 +600,7 @@ public:
 	 */
 	std::string toString() const {
 		std::stringstream result;
-		result << "bit vector(" << this->getNumberOfSetBits() << ") [";
+		result << "bit vector(" << this->getNumberOfSetBits() << "/" << bitCount << ") [";
 		for (auto index : *this) {
 			result << index << " ";
 		}
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 5e2019142..4ee005881 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -713,14 +713,18 @@ public:
 		SparseMatrix result(constraint.getNumberOfSetBits());
 		result.initialize(subNonZeroEntries);
 
-		// Create a temporary array that stores for each index whose bit is set
-		// to true the number of bits that were set before that particular index.
-		uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[colCount];
+		// Create a temporary vecotr that stores for each index whose bit is set
+        // to true the number of bits that were set before that particular index.
+        std::vector<uint_fast64_t> bitsSetBeforeIndex;
+        bitsSetBeforeIndex.reserve(colCount);
+        
+        // Compute the information to fill this vector.
 		uint_fast64_t lastIndex = 0;
 		uint_fast64_t currentNumberOfSetBits = 0;
 		for (auto index : constraint) {
 			while (lastIndex <= index) {
-				bitsSetBeforeIndex[lastIndex++] = currentNumberOfSetBits;
+				bitsSetBeforeIndex.push_back(currentNumberOfSetBits);
+                ++lastIndex;
 			}
 			++currentNumberOfSetBits;
 		}
@@ -737,9 +741,6 @@ public:
 			++rowCount;
 		}
 
-		// Dispose of the temporary array.
-		delete[] bitsSetBeforeIndex;
-
 		// Finalize sub-matrix and return result.
 		result.finalize();
 		LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
@@ -758,8 +759,7 @@ public:
 	SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> const& rowGroupIndices) const {
 		LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size).");
 
-		// First, we need to determine the number of non-zero entries and the number of rows of the
-		// sub-matrix.
+		// First, we need to determine the number of non-zero entries and the number of rows of the sub-matrix.
 		uint_fast64_t subNonZeroEntries = 0;
 		uint_fast64_t subRowCount = 0;
 		for (auto index : rowGroupConstraint) {
@@ -779,14 +779,18 @@ public:
 		SparseMatrix result(subRowCount, rowGroupConstraint.getNumberOfSetBits());
 		result.initialize(subNonZeroEntries);
 
-		// Create a temporary array that stores for each index whose bit is set
+		// Create a temporary vector that stores for each index whose bit is set
 		// to true the number of bits that were set before that particular index.
-		uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[colCount];
+        std::vector<uint_fast64_t> bitsSetBeforeIndex;
+        bitsSetBeforeIndex.reserve(colCount);
+        
+        // Compute the information to fill this vector.
 		uint_fast64_t lastIndex = 0;
 		uint_fast64_t currentNumberOfSetBits = 0;
 		for (auto index : rowGroupConstraint) {
 			while (lastIndex <= index) {
-				bitsSetBeforeIndex[lastIndex++] = currentNumberOfSetBits;
+				bitsSetBeforeIndex.push_back(currentNumberOfSetBits);
+                ++lastIndex;
 			}
 			++currentNumberOfSetBits;
 		}
@@ -804,9 +808,6 @@ public:
 			}
 		}
 
-		// Dispose of the temporary array.
-		delete[] bitsSetBeforeIndex;
-
 		// Finalize sub-matrix and return result.
 		result.finalize();
 		LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
diff --git a/src/utility/graph.h b/src/utility/graph.h
index dc8980b80..8d9f0f024 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -31,10 +31,12 @@ namespace graph {
      * @param model The model whose graph structure to search.
 	 * @param phiStates A bit vector of all states satisfying phi.
 	 * @param psiStates A bit vector of all states satisfying psi.
+     * @param useStepBound A flag that indicates whether or not to use the given number of maximal steps for the search.
+     * @param maximalSteps The maximal number of steps to reach the psi states. 
 	 * @return A bit vector with all indices of states that have a probability greater than 0.
 	 */
 	template <typename T>
-	storm::storage::BitVector performProbGreater0(storm::models::AbstractDeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+	storm::storage::BitVector performProbGreater0(storm::models::AbstractDeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) {
         // Prepare the resulting bit vector.
         storm::storage::BitVector statesWithProbabilityGreater0(model.getNumberOfStates());
         
@@ -48,16 +50,43 @@ namespace graph {
 		std::vector<uint_fast64_t> stack;
 		stack.reserve(model.getNumberOfStates());
 		psiStates.addSetIndicesToVector(stack);
+        
+        // Initialize the stack for the step bound, if the number of steps is bounded.
+        std::vector<uint_fast64_t> stepStack;
+        std::vector<uint_fast64_t> remainingSteps;
+        if (useStepBound) {
+            stepStack.reserve(model.getNumberOfStates());
+            stepStack.insert(stepStack.begin(), psiStates.getNumberOfSetBits(), maximalSteps);
+            remainingSteps.resize(model.getNumberOfStates());
+            for (auto state : psiStates) {
+                remainingSteps[state] = maximalSteps;
+            }
+        }
 
 		// Perform the actual BFS.
-		while(!stack.empty()) {
-			uint_fast64_t currentState = stack.back();
+        uint_fast64_t currentState, currentStepBound;
+		while (!stack.empty()) {
+            currentState = stack.back();
 			stack.pop_back();
+            
+            if (useStepBound) {
+                currentStepBound = stepStack.back();
+                stepStack.pop_back();
+            }
 
-			for(auto it = backwardTransitions.constColumnIteratorBegin(currentState); it != backwardTransitions.constColumnIteratorEnd(currentState); ++it) {
-				if (phiStates.get(*it) && !statesWithProbabilityGreater0.get(*it)) {
-					statesWithProbabilityGreater0.set(*it, true);
-					stack.push_back(*it);
+			for (auto it = backwardTransitions.constColumnIteratorBegin(currentState); it != backwardTransitions.constColumnIteratorEnd(currentState); ++it) {
+				if (phiStates.get(*it) && (!statesWithProbabilityGreater0.get(*it) || (useStepBound && remainingSteps[*it] < currentStepBound - 1))) {
+                    // If we don't have a number of maximal steps to take, just add the state to the stack.
+                    if (!useStepBound) {
+                        statesWithProbabilityGreater0.set(*it, true);
+                        stack.push_back(*it);
+                    } else if (currentStepBound > 0) {
+                        // If there is at least one more step to go, we need to push the state and the new number of steps.
+                        remainingSteps[*it] = currentStepBound - 1;
+                        statesWithProbabilityGreater0.set(*it, true);
+                        stack.push_back(*it);
+                        stepStack.push_back(currentStepBound - 1);
+                    }
 				}
 			}
 		}
@@ -127,45 +156,97 @@ namespace graph {
 	}
 
     /*!
-	 * Computes the sets of states that have probability 0 of satisfying phi until psi under all
-     * possible resolutions of non-determinism in a non-deterministic model. Stated differently,
-     * this means that these states have probability 0 of satisfying phi until psi even if the
-     * scheduler tries to maximize this probability.
+	 * Computes the sets of states that have probability greater 0 of satisfying phi until psi under at least
+     * one possible resolution of non-determinism in a non-deterministic model. Stated differently,
+     * this means that these states have a probability greater 0 of satisfying phi until psi if the
+     * scheduler tries to minimize this probability.
      *
 	 * @param model The model whose graph structure to search.
      * @param backwardTransitions The reversed transition relation of the model.
 	 * @param phiStates The set of all states satisfying phi.
 	 * @param psiStates The set of all states satisfying psi.
+     * @param useStepBound A flag that indicates whether or not to use the given number of maximal steps for the search.
+     * @param maximalSteps The maximal number of steps to reach the psi states.
 	 * @return A bit vector that represents all states with probability 0.
 	 */
 	template <typename T>
-	storm::storage::BitVector performProb0A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
-        // Prepare the resulting bit vector.
-        storm::storage::BitVector statesWithProbability0(model.getNumberOfStates());
-
+	storm::storage::BitVector performProbGreater0E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) {
+        // Prepare resulting bit vector.
+        storm::storage::BitVector statesWithProbabilityGreater0(model.getNumberOfStates());
+        
+		// Get some temporaries for convenience.
+		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
+		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+        
 		// Add all psi states as the already satisfy the condition.
-		statesWithProbability0 |= psiStates;
-
+		statesWithProbabilityGreater0 |= psiStates;
+        
 		// Initialize the stack used for the BFS with the states
 		std::vector<uint_fast64_t> stack;
 		stack.reserve(model.getNumberOfStates());
 		psiStates.addSetIndicesToVector(stack);
-
+        
+        // Initialize the stack for the step bound, if the number of steps is bounded.
+        std::vector<uint_fast64_t> stepStack;
+        std::vector<uint_fast64_t> remainingSteps;
+        if (useStepBound) {
+            stepStack.reserve(model.getNumberOfStates());
+            stepStack.insert(stepStack.begin(), psiStates.getNumberOfSetBits(), maximalSteps);
+            remainingSteps.resize(model.getNumberOfStates());
+            for (auto state : psiStates) {
+                remainingSteps[state] = maximalSteps;
+            }
+        }
+        
 		// Perform the actual BFS.
-		while(!stack.empty()) {
-			uint_fast64_t currentState = stack.back();
+        uint_fast64_t currentState, currentStepBound;
+		while (!stack.empty()) {
+			currentState = stack.back();
 			stack.pop_back();
-
-			for(auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
-				if (phiStates.get(*it) && !statesWithProbability0.get(*it)) {
-					statesWithProbability0.set(*it, true);
-					stack.push_back(*it);
-				}
+            
+            if (useStepBound) {
+                currentStepBound = stepStack.back();
+                stepStack.pop_back();
+            }
+            
+			for (auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
+                if (phiStates.get(*it) && (!statesWithProbabilityGreater0.get(*it) || (useStepBound && remainingSteps[*it] < currentStepBound - 1))) {
+                    // If we don't have a bound on the number of steps to take, just add the state to the stack.
+                    if (!useStepBound) {
+                        statesWithProbabilityGreater0.set(*it, true);
+                        stack.push_back(*it);
+                    } else if (currentStepBound > 0) {
+                        // If there is at least one more step to go, we need to push the state and the new number of steps.
+                        remainingSteps[*it] = currentStepBound - 1;
+                        statesWithProbabilityGreater0.set(*it, true);
+                        stack.push_back(*it);
+                        stepStack.push_back(currentStepBound - 1);
+                    }
+                }
 			}
 		}
-
-        // Finally, invert the computed set of states and return result.
-		statesWithProbability0.complement();
+        
+        return statesWithProbabilityGreater0;
+	}
+    
+    /*!
+	 * Computes the sets of states that have probability 0 of satisfying phi until psi under all
+     * possible resolutions of non-determinism in a non-deterministic model. Stated differently,
+     * this means that these states have probability 0 of satisfying phi until psi even if the
+     * scheduler tries to maximize this probability.
+     *
+	 * @param model The model whose graph structure to search.
+     * @param backwardTransitions The reversed transition relation of the model.
+	 * @param phiStates The set of all states satisfying phi.
+	 * @param psiStates The set of all states satisfying psi.
+     * @param useStepBound A flag that indicates whether or not to use the given number of maximal steps for the search.
+     * @param maximalSteps The maximal number of steps to reach the psi states. 
+	 * @return A bit vector that represents all states with probability 0.
+	 */
+	template <typename T>
+	storm::storage::BitVector performProb0A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+        storm::storage::BitVector statesWithProbability0 = performProbGreater0E(model, backwardTransitions, phiStates, psiStates);
+        statesWithProbability0.complement();
         return statesWithProbability0;
 	}
 
@@ -194,13 +275,14 @@ namespace graph {
 
         // Perform the loop as long as the set of states gets larger.
 		bool done = false;
+        uint_fast64_t currentState;
 		while (!done) {
 			stack.clear();
 			storm::storage::BitVector nextStates(psiStates);
 			psiStates.addSetIndicesToVector(stack);
 
 			while (!stack.empty()) {
-				uint_fast64_t currentState = stack.back();
+				currentState = stack.back();
 				stack.pop_back();
 
 				for(auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
@@ -263,69 +345,116 @@ namespace graph {
 	}
 
     /*!
-	 * Computes the sets of states that have probability 0 of satisfying phi until psi under at least
-     * one possible resolution of non-determinism in a non-deterministic model. Stated differently,
-     * this means that these states have probability 0 of satisfying phi until psi if the
-     * scheduler tries to minimize this probability.
+	 * Computes the sets of states that have probability greater 0 of satisfying phi until psi under any
+     * possible resolution of non-determinism in a non-deterministic model. Stated differently,
+     * this means that these states have a probability greater 0 of satisfying phi until psi if the
+     * scheduler tries to maximize this probability.
      *
 	 * @param model The model whose graph structure to search.
      * @param backwardTransitions The reversed transition relation of the model.
 	 * @param phiStates The set of all states satisfying phi.
 	 * @param psiStates The set of all states satisfying psi.
+     * @param useStepBound A flag that indicates whether or not to use the given number of maximal steps for the search.
+     * @param maximalSteps The maximal number of steps to reach the psi states.
 	 * @return A bit vector that represents all states with probability 0.
 	 */
 	template <typename T>
-	storm::storage::BitVector performProb0E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+	storm::storage::BitVector performProbGreater0A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) {
         // Prepare resulting bit vector.
-        storm::storage::BitVector statesWithProbability0(model.getNumberOfStates());
+        storm::storage::BitVector statesWithProbabilityGreater0(model.getNumberOfStates());
         
 		// Get some temporaries for convenience.
 		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
 		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
-
+        
 		// Add all psi states as the already satisfy the condition.
-		statesWithProbability0 |= psiStates;
-
-		// Initialize the stack used for the DFS with the states
+		statesWithProbabilityGreater0 |= psiStates;
+        
+		// Initialize the stack used for the BFS with the states
 		std::vector<uint_fast64_t> stack;
 		stack.reserve(model.getNumberOfStates());
 		psiStates.addSetIndicesToVector(stack);
-
-		// Perform the actual DFS.
+        
+        // Initialize the stack for the step bound, if the number of steps is bounded.
+        std::vector<uint_fast64_t> stepStack;
+        std::vector<uint_fast64_t> remainingSteps;
+        if (useStepBound) {
+            stepStack.reserve(model.getNumberOfStates());
+            stepStack.insert(stepStack.begin(), psiStates.getNumberOfSetBits(), maximalSteps);
+            remainingSteps.resize(model.getNumberOfStates());
+            for (auto state : psiStates) {
+                remainingSteps[state] = maximalSteps;
+            }
+        }
+        
+		// Perform the actual BFS.
+        uint_fast64_t currentState, currentStepBound;
 		while(!stack.empty()) {
-			uint_fast64_t currentState = stack.back();
+			currentState = stack.back();
 			stack.pop_back();
-
+            
+            if (useStepBound) {
+                currentStepBound = stepStack.back();
+                stepStack.pop_back();
+            }
+            
 			for(auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
-				if (phiStates.get(*it) && !statesWithProbability0.get(*it)) {
-					// Check whether the predecessor has at least one successor in the current state
-					// set for every nondeterministic choice.
-					bool addToStatesWithProbability0 = true;
-					for (auto rowIt = nondeterministicChoiceIndices->begin() + *it; rowIt != nondeterministicChoiceIndices->begin() + *it + 1; ++rowIt) {
-						bool hasAtLeastOneSuccessorWithProbabilityGreater0 = false;
-						for (auto colIt = transitionMatrix->constColumnIteratorBegin(*rowIt); colIt != transitionMatrix->constColumnIteratorEnd(*rowIt); ++colIt) {
-							if (statesWithProbability0.get(*colIt)) {
-								hasAtLeastOneSuccessorWithProbabilityGreater0 = true;
-								break;
-							}
-						}
-						if (!hasAtLeastOneSuccessorWithProbabilityGreater0) {
-							addToStatesWithProbability0 = false;
-							break;
-						}
-					}
-
-					// If we need to add the state, then actually add it and perform further search
-					// from the state.
-					if (addToStatesWithProbability0) {
-						statesWithProbability0.set(*it, true);
-						stack.push_back(*it);
-					}
-				}
+                if (phiStates.get(*it) && (!statesWithProbabilityGreater0.get(*it) || (useStepBound && remainingSteps[*it] < currentStepBound - 1))) {
+                    // Check whether the predecessor has at least one successor in the current state set for every
+                    // nondeterministic choice.
+                    bool addToStatesWithProbabilityGreater0 = true;
+                    for (auto rowIt = nondeterministicChoiceIndices->begin() + *it; rowIt != nondeterministicChoiceIndices->begin() + *it + 1; ++rowIt) {
+                        bool hasAtLeastOneSuccessorWithProbabilityGreater0 = false;
+                        for (auto colIt = transitionMatrix->constColumnIteratorBegin(*rowIt); colIt != transitionMatrix->constColumnIteratorEnd(*rowIt); ++colIt) {
+                            if (statesWithProbabilityGreater0.get(*colIt)) {
+                                hasAtLeastOneSuccessorWithProbabilityGreater0 = true;
+                                break;
+                            }
+                        }
+                        
+                        if (!hasAtLeastOneSuccessorWithProbabilityGreater0) {
+                            addToStatesWithProbabilityGreater0 = false;
+                            break;
+                        }
+                    }
+                    
+                    // If we need to add the state, then actually add it and perform further search from the state.
+                    if (addToStatesWithProbabilityGreater0) {
+                        // If we don't have a bound on the number of steps to take, just add the state to the stack.
+                        if (!useStepBound) {
+                            statesWithProbabilityGreater0.set(*it, true);
+                            stack.push_back(*it);
+                        } else if (currentStepBound > 0) {
+                            // If there is at least one more step to go, we need to push the state and the new number of steps.
+                            remainingSteps[*it] = currentStepBound - 1;
+                            statesWithProbabilityGreater0.set(*it, true);
+                            stack.push_back(*it);
+                            stepStack.push_back(currentStepBound - 1);
+                        }
+                    }
+                }
 			}
 		}
-
-		statesWithProbability0.complement();
+        
+        return statesWithProbabilityGreater0;
+	}
+    
+    /*!
+	 * Computes the sets of states that have probability 0 of satisfying phi until psi under at least
+     * one possible resolution of non-determinism in a non-deterministic model. Stated differently,
+     * this means that these states have probability 0 of satisfying phi until psi if the
+     * scheduler tries to minimize this probability.
+     *
+	 * @param model The model whose graph structure to search.
+     * @param backwardTransitions The reversed transition relation of the model.
+	 * @param phiStates The set of all states satisfying phi.
+	 * @param psiStates The set of all states satisfying psi.
+	 * @return A bit vector that represents all states with probability 0.
+	 */
+	template <typename T>
+	storm::storage::BitVector performProb0E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+        storm::storage::BitVector statesWithProbability0 = performProbGreater0A(model, backwardTransitions, phiStates, psiStates);
+        statesWithProbability0.complement();
         return statesWithProbability0;
 	}
 
@@ -347,20 +476,21 @@ namespace graph {
 		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
 		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
 
+        // Initialize the environment for the iterative algorithm.
 		storm::storage::BitVector currentStates(model.getNumberOfStates(), true);
-
 		std::vector<uint_fast64_t> stack;
 		stack.reserve(model.getNumberOfStates());
 
         // Perform the loop as long as the set of states gets smaller.
 		bool done = false;
+        uint_fast64_t currentState;
 		while (!done) {
 			stack.clear();
 			storm::storage::BitVector nextStates(psiStates);
 			psiStates.addSetIndicesToVector(stack);
 
 			while (!stack.empty()) {
-				uint_fast64_t currentState = stack.back();
+				currentState = stack.back();
 				stack.pop_back();
 
 				for(auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
index 6efb4653a..0d12e3a94 100644
--- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
@@ -209,7 +209,7 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
 	storm::property::prctl::BoundedEventually<double>* boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, true);
+	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, false);
 
 	result = mc.checkNoBoundOperator(*probFormula);
 
@@ -218,11 +218,11 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 	ASSERT_LT(std::abs((*result)[0] - 0.0625), s->get<double>("precision"));
 
 	delete probFormula;
-	delete result;
+	delete result; 
 
 	apFormula = new storm::property::prctl::Ap<double>("elected");
 	boundedEventuallyFormula = new storm::property::prctl::BoundedEventually<double>(apFormula, 25);
-	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, false);
+	probFormula = new storm::property::prctl::ProbabilisticNoBoundOperator<double>(boundedEventuallyFormula, true);
 
 	result = mc.checkNoBoundOperator(*probFormula);