diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index bd2c99760..de2459647 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -134,14 +134,23 @@ namespace storm {
             void SparseMarkovAutomatonCslHelper::computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory) {
                 STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded reachability probabilities is unsupported for this value type.");
             }
-            
-            template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
-            void SparseMarkovAutomatonCslHelper::printTransitions(std::vector<std::vector<ValueType>> relReachability, std::vector<ValueType> const& exitRateVector, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector const& cycleStates, storm::storage::BitVector const& cycleGoalStates, std::vector<std::vector<ValueType>>& vd, std::vector<std::vector<ValueType>>& vu, std::vector<std::vector<ValueType>>& wu){
+
+            template<typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type= 0>
+            void SparseMarkovAutomatonCslHelper::printTransitions(
+                    storm::storage::SparseMatrix<ValueType> const &fullTransitionMatrix,
+                    std::vector<ValueType> const &exitRateVector, storm::storage::BitVector const &markovianStates,
+                    storm::storage::BitVector const &psiStates, std::vector<std::vector<ValueType>> relReachability,
+                    const storage::BitVector &cycleStates, const storage::BitVector &cycleGoalStates,
+                    std::vector<std::vector<std::vector<ValueType>>> &unifVectors) {
+
+
+
                 std::ofstream logfile("U+logfile.txt", std::ios::app);
-                
+
                 auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices();
                 auto numberOfStates = fullTransitionMatrix.getRowGroupCount();            
 
+                //Transition Matrix
                 logfile << "number of states = num of row group count " << numberOfStates << "\n";
                 for (uint_fast64_t i = 0; i < fullTransitionMatrix.getRowGroupCount(); i++) {
                     logfile << " from node  " << i << " ";
@@ -161,8 +170,8 @@ namespace storm {
                 logfile << "\n";
 
                 logfile << "probStates\tmarkovianStates\tgoalStates\tcycleStates\tcycleGoalStates\n";
-                for (int i =0 ; i< markovianStates.size() ; i++){
-                    logfile << (~markovianStates)[i] << "\t\t" << markovianStates[i] << "\t\t" << psiStates[i] << "\t\t" << cycleStates[i] << "\t\t" << cycleGoalStates[i] << "\n";
+                    for (int i =0 ; i< markovianStates.size() ; i++){
+                        logfile << (~markovianStates)[i] << "\t\t" << markovianStates[i] << "\t\t" << psiStates[i] << "\t\t" << cycleStates[i] << "\t\t" << cycleGoalStates[i] << "\n";
                 }
 
                 for (int i =0; i<relReachability.size(); i++){
@@ -174,25 +183,25 @@ namespace storm {
 
 
                 logfile << "vd: \n";
-                for (uint_fast64_t i =0 ; i<vd.size(); i++){
-                    for(uint_fast64_t j=0; j<fullTransitionMatrix.getRowGroupCount(); j++){
-                        logfile << vd[i][j] << "\t" ;
+                for (uint64_t i =0 ; i<unifVectors[0].size(); i++){
+                    for(uint64_t j=0; j<unifVectors[0][i].size(); j++){
+                        logfile << unifVectors[0][i][j] << "\t" ;
                     }
                     logfile << "\n";
                 }
 
                 logfile << "\nvu:\n";
-                for (uint_fast64_t i =0 ; i<vu.size(); i++){
-                    for(uint_fast64_t j=0; j<fullTransitionMatrix.getRowGroupCount(); j++){
-                        logfile << vu[i][j] << "\t" ;
+                for (uint64_t i =0 ; i<unifVectors[1].size(); i++){
+                    for(uint64_t j=0; j<unifVectors[1][i].size(); j++){
+                        logfile << unifVectors[1][i][j] << "\t" ;
                     }
                     logfile << "\n";
                 }
 
                 logfile << "\nwu\n";
-                for (uint_fast64_t i =0 ; i<wu.size(); i++){
-                    for(uint_fast64_t j=0; j<fullTransitionMatrix.getRowGroupCount(); j++){
-                        logfile << wu[i][j] << "\t" ;
+                for (uint64_t i =0 ; i<unifVectors[2].size(); i++){
+                    for(uint64_t j=0; j<unifVectors[2][i].size(); j++){
+                        logfile << unifVectors[2][i][j] << "\t" ;
                     }
                     logfile << "\n";
                 }
@@ -212,256 +221,137 @@ namespace storm {
                 return res;
             }
 
+
+
             template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
-            void SparseMarkovAutomatonCslHelper::calculateVu(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector<std::vector<ValueType>>& vu, std::vector<std::vector<ValueType>>& wu, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver){
-                if (vu[k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation.
-                uint64_t N = vu.size()-1;
+            void SparseMarkovAutomatonCslHelper::calculateVu(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, uint64_t const kind, ValueType lambda, uint64_t probSize, std::vector<std::vector<std::vector<ValueType>>>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver){
+                if (unifVectors[1][k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation.
+                uint64_t N = unifVectors[1].size()-1;
                 auto rowGroupIndices = fullTransitionMatrix.getRowGroupIndices();
 
                 ValueType res =0;
                 for (uint64_t i = k ; i < N ; i++ ){
-                    if (wu[N-1-(i-k)][node]==-1){
-                        calculateWu(relativeReachability, dir, (N-1-(i-k)),node,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, solver);
+                    if (unifVectors[2][N-1-(i-k)][node]==-1){
+                       calculateUnifPlusVector(N-1-(i-k),node,2,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix, markovianStates,psiStates,solver);
+                               //old:  relativeReachability, dir, (N-1-(i-k)),node,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, solver);
                     }
-                    res+=poisson(lambda, i)*wu[N-1-(i-k)][node];
+                    res+=poisson(lambda, i)*unifVectors[2][N-1-(i-k)][node];
                 }
-                vu[k][node]=res;
+                unifVectors[1][k][node]=res;
             }
 
-            template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
-            void SparseMarkovAutomatonCslHelper::calculateWu(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector<std::vector<ValueType>>& wu, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver){
-                if (wu[k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation.
-
-                std::ofstream logfile("U+logfile.txt", std::ios::app);
-
-                auto probabilisticStates = ~markovianStates;
 
-                uint64_t N = wu.size()-1;
-                auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices( );
 
-                ValueType res;
-                if (k==N){
-                    wu[k][node]=0;
-                    return;
-                }
 
-                if (psiStates[node]){
-                    wu[k][node]=1;
-                    return;
-                }
+            template<typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type= 0>
+            void SparseMarkovAutomatonCslHelper::calculateUnifPlusVector(uint64_t k, uint64_t node, uint64_t const kind, ValueType lambda, uint64_t probSize,
+                                                                         std::vector<std::vector<ValueType>> const &relativeReachability,
+                                                                         OptimizationDirection dir,
+                                                                         std::vector<std::vector<std::vector<ValueType>>> &unifVectors,
+                                                                         storm::storage::SparseMatrix<ValueType> const &fullTransitionMatrix,
+                                                                         storm::storage::BitVector const &markovianStates,
+                                                                         storm::storage::BitVector const &psiStates,
+                                                                         std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const &solver) {
 
-                if (markovianStates[node]){
-                    res = 0;
-                    auto line = fullTransitionMatrix.getRow(rowGroupIndices[node]);
-                    for (auto &element : line){
-                        uint64_t to = element.getColumn();
-                        if (wu[k+1][to]==-1){
-                            calculateWu(relativeReachability, dir, k+1,to,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, solver);
-                        }
-                        res+=element.getValue()*wu[k+1][to];
-                    }
-                wu[k][node]=res;
-                } else {
-                    // applying simpleMDP stuff
-                    auto numberOfStates=relativeReachability[0].size();
-                    typename storm::storage::SparseMatrix<ValueType> probMatrix = fullTransitionMatrix.getSubmatrix(true, ~markovianStates , ~markovianStates, true);
-                    auto probSize = probMatrix.getRowGroupCount();
-                    std::vector<ValueType> b(probMatrix.getRowCount(), 0), x(probSize,0);
-
-                    //calculate b
-                    uint64_t  lineCounter=0;
-                    for (int i =0; i<numberOfStates; i++){
-                        if (markovianStates[i]){
-                            continue;
-                        }
-                        auto rowStart = rowGroupIndices[i];
-                        auto rowEnd = rowGroupIndices[i+1];
-                        for (auto j = rowStart; j<rowEnd; j++){
-                            ValueType res=0;
-                            for (auto& element:fullTransitionMatrix.getRow(j)){
-                                auto to = element.getColumn();
-                                if (probabilisticStates[to]){
-                                    continue;
-                                }
-                                if (wu[k][to]==-1){
-                                    calculateWu(relativeReachability, dir, k, to, lambda, wu, fullTransitionMatrix, markovianStates, psiStates, solver);
-                                }
-                                res = res + relativeReachability[j][to]*wu[k][to];
-                            }
-                            b[lineCounter]=res;
-                            lineCounter++;
-                        }
-                    }
-                    logfile << "\nwu reach vector b is \n";
-                    for (uint64_t i =0 ; i<probSize; i++){
-                        logfile << b[i] << "\n";
-                    }
-
-                    solver->solveEquations(dir, x, b);
-
-                    for (uint64_t i =0 ; i<probSize ; i++){
-                        auto trueI = transformIndice(~markovianStates,i);
-                        wu[k][trueI]=x[i];
-                    }
-
-                    /* unifPlus
-                    res = -1;
-                    uint64_t rowStart = rowGroupIndices[node];
-                    uint64_t rowEnd = rowGroupIndices[node+1];
-                    for (uint64_t i = rowStart; i< rowEnd; i++){
-                        auto line = fullTransitionMatrix.getRow(i);
-                        ValueType between = 0;
-                        for (auto& element: line){
-                            uint64_t to = element.getColumn();
-                            if (to==node){
-                                continue;
-                            }
-                            if (wu[k][to]==-1){
-                                calculateWu(dir, k,to,lambda,wu,fullTransitionMatrix,markovianStates,psiStates);
-                            }
-                            between+=element.getValue()*wu[k][to];
-                        }
-                        if (maximize(dir)){
-                            res = std::max(res,between);
-                        } else {
-                            if (res!=-1){
-                                res = std::min(res,between);
-                            } else {
-                                res = between;
-                            }
-                        }
-                    } */
-                } // end no goal-prob state
-            }
-
-            template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
-            void SparseMarkovAutomatonCslHelper::calculateVd(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector<std::vector<ValueType>>& vd, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver){
-               
+                if (unifVectors[kind][k][node]!=-1){return;}
                 std::ofstream logfile("U+logfile.txt", std::ios::app);
+                logfile << "calculating vector " << kind   <<  " for k = " << k << " node "<< node << " \t";
 
                 auto probabilisticStates = ~markovianStates;
-
-                if (vd[k][node]!=-1){return;} //dynamic programming. avoiding multiple calculation.
-                logfile << "calculating vd for k = " << k << " node "<< node << " \t";
-                uint64_t N = vd.size()-1;
+                auto numberOfStates=fullTransitionMatrix.getRowGroupCount();
+                uint64_t N = unifVectors[kind].size()-1;
                 auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices();
-
                 ValueType res;
+
+                // First Case, k==N, independent from kind of state
                 if (k==N){
                     logfile << "k == N! res = 0\n";
-                    vd[k][node]=0;
+                    unifVectors[kind][k][node]=0;
                     return;
                 }
-                
+
                 //goal state
                 if (psiStates[node]){
-                    res = storm::utility::zero<ValueType>();
-                    for (uint64_t i = k ; i<N ; i++){
-                        ValueType between = poisson(lambda,i);
-                        res+=between;
+                    if (kind==0){
+                        // Vd
+                        res = storm::utility::zero<ValueType>();
+                        for (uint64_t i = k ; i<N ; i++){
+                            ValueType between = poisson(lambda,i);
+                            res+=between;
+                        }
+                        unifVectors[kind][k][node]=res;
+                    } else {
+                        // WU
+                        unifVectors[kind][k][node]=1;
                     }
-                    vd[k][node]=res;
                     logfile << "goal state node " << node << " res = " << res << "\n";
                     return;
+
                 }
 
-                // no-goal markovian state
+                //markovian non-goal State
                 if (markovianStates[node]){
                     logfile << "markovian state: ";
                     res = 0;
                     auto line = fullTransitionMatrix.getRow(rowGroupIndices[node]);
                     for (auto &element : line){
                         uint64_t to = element.getColumn();
-                        if (vd[k+1][to]==-1){
-                            calculateVd(relativeReachability, dir,k+1,to,lambda,vd, fullTransitionMatrix, markovianStates,psiStates, solver);
+                        if (unifVectors[kind][k+1][to]==-1){
+                            calculateUnifPlusVector(k+1,to,kind,lambda,probSize,relativeReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver);
                         }
-                        res+=element.getValue()*vd[k+1][to];
+                        res+=element.getValue()*unifVectors[kind][k+1][to];
                     }
-
-                    vd[k][node]=res;
+                    unifVectors[kind][k][node]=res;
                     logfile << " res = " << res << "\n";
-                } else { //no-goal prob state
-                    // applying simpleMDP stuff
-                    logfile << "probState, applying MDP stuff ";
-                    auto numberOfStates=relativeReachability[0].size();
-                    typename storm::storage::SparseMatrix<ValueType> probMatrix = fullTransitionMatrix.getSubmatrix(true, ~markovianStates , ~markovianStates, true);
-                    auto probSize = probMatrix.getRowGroupCount();
-                    std::vector<ValueType> b(probMatrix.getRowCount(), 0), x(probSize,0);
+                    return;
+                }
 
+                //probabilistic non-goal State
+                if (probabilisticStates[node]){
+                    logfile << "probabilistic state: ";
+                    std::vector<ValueType> b(probSize, 0), x(probabilisticStates.getNumberOfSetBits(),0);
                     //calculate b
                     uint64_t  lineCounter=0;
-                    for (int i =0; i<numberOfStates; i++){
-                        if (markovianStates[i]){
+                    for (int i =0; i<numberOfStates; i++) {
+                        if (markovianStates[i]) {
                             continue;
                         }
                         auto rowStart = rowGroupIndices[i];
-                        auto rowEnd = rowGroupIndices[i+1];
-                        for (auto j = rowStart; j<rowEnd; j++){
-                            ValueType res=0;
-                            for (auto& element:fullTransitionMatrix.getRow(j)){
+                        auto rowEnd = rowGroupIndices[i + 1];
+                        for (auto j = rowStart; j < rowEnd; j++) {
+                            res = 0;
+                            for (auto &element:fullTransitionMatrix.getRow(j)) {
                                 auto to = element.getColumn();
-                                if (probabilisticStates[to]){
+                                if (probabilisticStates[to]) {
                                     continue;
                                 }
-                                if (vd[k][to]==-1){
-                                    calculateVd(relativeReachability, dir, k, to, lambda, vd, fullTransitionMatrix, markovianStates, psiStates, solver);
+                                if (unifVectors[kind][k][to] == -1) {
+                                    calculateUnifPlusVector(k, to, kind, lambda, probSize, relativeReachability, dir,
+                                                            unifVectors, fullTransitionMatrix, markovianStates,
+                                                            psiStates, solver);
                                 }
-                                res = res + relativeReachability[j][to]*vd[k][to];
+                                res = res + relativeReachability[j][to] * unifVectors[kind][k][to];
                             }
-                            b[lineCounter]=res;
+                            b[lineCounter] = res;
                             lineCounter++;
                         }
                     }
 
-                    logfile << "\nvd reach vector b is \n";
-                    for (uint64_t i =0 ; i<probMatrix.getRowCount(); i++){
-                        logfile << b[i] << "\n";
-                    }
-
                     solver->solveEquations(dir, x, b);
-                    logfile << "vd goal vector x is \n";
-                    for (int i =0 ; i<probMatrix.getRowCount(); i++){
-                        logfile << x[i] << "\n";
-                    }
+
 
 
                     for (uint64_t i =0 ; i<probSize ; i++){
-                        auto trueI = transformIndice(~markovianStates,i);
-                        vd[k][trueI]=x[i];
-                    }
-                    /* unif Plus
-                    logfile << "prob state: ";
-                    res = -1;
-                    uint64_t rowStart = rowGroupIndices[node];
-                    uint64_t rowEnd = rowGroupIndices[node+1];
-                    for (uint64_t i = rowStart; i< rowEnd; i++){
-                        auto line = fullTransitionMatrix.getRow(i);
-                        ValueType between = 0;
-                        for (auto& element: line){
-                            uint64_t to = element.getColumn();
-                            if (to==node){
-                                logfile << "ignoring self loops for now";
-                                continue;
-                            }
-                            if (vd[k][to]==-1){
-                                calculateVd(dir, k,to,lambda,vd, fullTransitionMatrix, markovianStates,psiStates);
-                            }
-                            between+=element.getValue()*vd[k][to];
-                        }
-                        if (maximize(dir)){
-                                res = std::max(res, between);
-                        } else {
-                            if (res!=-1){
-                                res =std::min(res,between);
-                            } else {
-                                res = between;
-                            }
-                        }
+                        auto trueI = transformIndice(probabilisticStates,i);
+                        unifVectors[kind][k][trueI]=x[i];
                     }
-                    */
-                }
+
+                    logfile << " res = " << unifVectors[kind][k][node] << " but calculated more \n";
+
+                } //end probabilistic states
             }
 
+
             template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
             uint64_t SparseMarkovAutomatonCslHelper::trajans(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, uint64_t node, std::vector<uint64_t >& disc, std::vector<uint64_t >& finish, uint64_t* counter) {
                 auto const& rowGroupIndice = transitionMatrix.getRowGroupIndices();
@@ -585,24 +475,29 @@ namespace storm {
             std::vector<ValueType> SparseMarkovAutomatonCslHelper::unifPlus(OptimizationDirection dir, std::pair<double, double> const& boundsPair, std::vector<ValueType> const& exitRateVector, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory){
                 STORM_LOG_TRACE("Using UnifPlus to compute bounded until probabilities.");
 
-                std::ofstream logfile("U+logfile.txt", std::ios::app);
-                ValueType maxNorm = 0;
+                ValueType maxNorm = storm::utility::zero<ValueType>();
 
                 //bitvectors to identify different kind of states
                 storm::storage::BitVector allStates(markovianStates.size(), true);
-                auto probabilisticStates = ~markovianStates;
+                storm::storage::BitVector probabilisticStates = ~markovianStates;
 
 
                 //vectors to save calculation
                 std::vector<std::vector<ValueType>> vd,vu,wu;
+                std::vector<std::vector<std::vector<ValueType>>> unifVectors{};
 
 
                 //transition matrix with diagonal entries. The values can be changed during uniformisation
                 std::vector<ValueType> exitRate{exitRateVector};
                 typename storm::storage::SparseMatrix<ValueType> fullTransitionMatrix = transitionMatrix.getSubmatrix(true, allStates , allStates , true);
-                typename storm::storage::SparseMatrix<ValueType> probMatrix = fullTransitionMatrix.getSubmatrix(true, ~markovianStates , ~markovianStates, true);
+                typename storm::storage::SparseMatrix<ValueType> probMatrix{};
+                uint64_t  probSize =0;
+                if (probabilisticStates.getNumberOfSetBits()!=0){ //work around in case there are no prob states
+                    probMatrix = fullTransitionMatrix.getSubmatrix(true, probabilisticStates , probabilisticStates, true);
+                    probSize = probMatrix.getRowCount();
+                }
 
-                auto rowGroupIndices = fullTransitionMatrix.getRowGroupIndices();
+                auto& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices();
 
 
                 //(1) define horizon, epsilon, kappa , N, lambda,
@@ -616,13 +511,14 @@ namespace storm {
                 }
                 uint64_t N;
 
+
+
                 //calculate relative ReachabilityVectors
                 std::vector<ValueType> in(numberOfStates, 0);
                 std::vector<std::vector<ValueType>> relReachability(fullTransitionMatrix.getRowCount(),in);
 
-                printTransitions(relReachability, exitRateVector, fullTransitionMatrix, markovianStates, psiStates, psiStates, psiStates, vd,vu,wu); // TODO: delete when develepmont is finished
 
-                 //calculate relative reachability
+                //calculate relative reachability
 
                 for(uint64_t i=0 ; i<numberOfStates; i++){
                     if (markovianStates[i]){
@@ -645,15 +541,14 @@ namespace storm {
                 requirements.clearBounds();
                 STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver.");
 
-                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(probMatrix);
-
-                solver->setHasUniqueSolution();
-                solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
-                solver->setRequirementsChecked();
-                solver->setCachingEnabled(true);
-
-
-                printTransitions(relReachability, exitRateVector, fullTransitionMatrix, markovianStates, psiStates, psiStates, psiStates, vd,vu,wu); // TODO: delete when develepmont is finished
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver;
+                if (probSize!=0){
+                    solver = minMaxLinearEquationSolverFactory.create(probMatrix);
+                    solver->setHasUniqueSolution();
+                    solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
+                    solver->setRequirementsChecked();
+                    solver->setCachingEnabled(true);
+                }
                     // while not close enough to precision:
                 do {
                     // (2) update parameter
@@ -692,29 +587,37 @@ namespace storm {
                     vd = std::vector<std::vector<ValueType>> (N + 1, init);
                     vu = std::vector<std::vector<ValueType>> (N + 1, init);
                     wu = std::vector<std::vector<ValueType>> (N + 1, init);
-                    
 
+                    unifVectors.clear();
+                    unifVectors.push_back(vd);
+                    unifVectors.push_back(vd);
+                    unifVectors.push_back(vd);
 
+
+                   printTransitions(fullTransitionMatrix,exitRate,markovianStates,psiStates,relReachability,psiStates, psiStates,unifVectors); //TODO remove
+
+                    //define 0=vd 1=vu 2=wu
                     // (5) calculate vectors and maxNorm
                     for (uint64_t i = 0; i < numberOfStates; i++) {
                         for (uint64_t k = N; k <= N; k--) {
-                                calculateVd(relReachability, dir, k, i, T*lambda, vd, fullTransitionMatrix, markovianStates, psiStates, solver);
-                                calculateWu(relReachability, dir, k, i, T*lambda, wu, fullTransitionMatrix, markovianStates, psiStates, solver);
-                                calculateVu(relReachability, dir, k, i, T*lambda, vu, wu, fullTransitionMatrix, markovianStates, psiStates, solver);
-                                //also use iteration to keep maxNorm of vd and vu up to date, so the loop-condition is easy to prove
-                                ValueType diff = std::abs(vd[k][i]-vu[k][i]);
+                                calculateUnifPlusVector(k,i,0,lambda,probSize,relReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver);
+                                calculateUnifPlusVector(k,i,2,lambda,probSize,relReachability,dir,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver);
+                                calculateVu(relReachability,dir,k,i,1,lambda,probSize,unifVectors,fullTransitionMatrix,markovianStates,psiStates,solver);
+                                //also use iteration to keep maxNorm of vd and vup to date, so the loop-condition is easy to prove
+                                ValueType diff = std::abs(unifVectors[0][k][i]-unifVectors[1][k][i]);
                                 maxNorm  = std::max(maxNorm, diff);
                             }
                     }
-
-                    printTransitions(relReachability, exitRateVector, fullTransitionMatrix, markovianStates, psiStates, psiStates, psiStates, vd,vu,wu); // TODO: delete when develepmont is finished
+                    std::cout << "\nTBU was " << unifVectors[0][0][0] << "\n";
+                    printTransitions(fullTransitionMatrix,exitRate,markovianStates,psiStates,relReachability,psiStates, psiStates,unifVectors); //TODO remove
 
 
-                    // (6) double lambda 
+                    // (6) double lambda
                     lambda=2*lambda;
 
                 } while (maxNorm>epsilon*(1-kappa));
-                return vd[0];
+
+                return unifVectors[0][0];
             }
 
             template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
@@ -1234,7 +1137,7 @@ namespace storm {
                 return v.front() * uniformizationRate;
             
             }
-            
+
             template std::vector<double> SparseMarkovAutomatonCslHelper::computeBoundedUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix<double> const& transitionMatrix, std::vector<double> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair<double, double> const& boundsPair, storm::solver::MinMaxLinearEquationSolverFactory<double> const& minMaxLinearEquationSolverFactory);
                 
             template std::vector<double> SparseMarkovAutomatonCslHelper::computeUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix<double> const& transitionMatrix, storm::storage::SparseMatrix<double> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory<double> const& minMaxLinearEquationSolverFactory);
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
index d378ded43..3d5b8d049 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
@@ -56,6 +56,8 @@ namespace storm {
                 static std::vector<ValueType> computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
                 
             private:
+                template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
+                static void calculateUnifPlusVector(uint64_t k, uint64_t node, uint64_t const kind, ValueType lambda, uint64_t probSize, std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, std::vector<std::vector<std::vector<ValueType>>>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver);
 
                 template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
                 static void deleteProbDiagonals(storm::storage::SparseMatrix<ValueType>& transitionMatrix, storm::storage::BitVector const& markovianStates);
@@ -93,27 +95,15 @@ namespace storm {
                 template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
                 static uint64_t trajans(storm::storage::SparseMatrix<ValueType> const& TransitionMatrix, uint64_t node, std::vector<uint64_t>& disc, std::vector<uint64_t>& finish, uint64_t * counter);
 
-                /*!
-                 * Computes vd vector according to UnifPlus
-                 *
-                 */
-                template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
-                static void calculateVd(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector<std::vector<ValueType>>& vd, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver);
-
-/*!
+                /*
                  * Computes vu vector according to UnifPlus
                  *
                  */
                 template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
-                static void calculateVu(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector<std::vector<ValueType>>& vu, std::vector<std::vector<ValueType>>& wu, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver);
+                static void calculateVu(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, uint64_t const kind, ValueType lambda, uint64_t probSize, std::vector<std::vector<std::vector<ValueType>>>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver);
+
 
 
-                /*!
-                 * Computes wu vector according to UnifPlus
-                 *
-                 */
-                template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
-                static void calculateWu(std::vector<std::vector<ValueType>> const& relativeReachability, OptimizationDirection dir, uint64_t k, uint64_t node, ValueType lambda, std::vector<std::vector<ValueType>>& wu, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates,  std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver);
 
                 /*!
                  * Prints the TransitionMatrix and the vectors vd, vu, wu to the logfile
@@ -122,7 +112,9 @@ namespace storm {
                  */
 
                 template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
-                static void printTransitions(std::vector<std::vector<ValueType>> relReachability, std::vector<ValueType> const& exitRateVector, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates,  storm::storage::BitVector const& cycleStates, storm::storage::BitVector const& cycleGoalStates, std::vector<std::vector<ValueType>>& vd, std::vector<std::vector<ValueType>>& vu, std::vector<std::vector<ValueType>>& wu);
+                static void printTransitions(storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, std::vector<ValueType> const& exitRateVector,  storm::storage::BitVector const& markovianStates,
+                                             storm::storage::BitVector const& psiStates,  std::vector<std::vector<ValueType>> relReachability,
+                                             storm::storage::BitVector const& cycleStates , storm::storage::BitVector const& cycleGoalStates ,std::vector<std::vector<std::vector<ValueType>>>& unifVectors);
 
                 template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type = 0>
                 static void computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);