diff --git a/ConSol/ConSol/Base/runSolver.cpp b/ConSol/ConSol/Base/runSolver.cpp
index 67c520b..14152a9 100644
--- a/ConSol/ConSol/Base/runSolver.cpp
+++ b/ConSol/ConSol/Base/runSolver.cpp
@@ -1,95 +1,110 @@
 //
 //  runSolver.cpp
 //  ConSol
 //
 //  Created by Fabian Moenkeberg on 30.06.17.
 //  Copyright © 2017 EPFL. All rights reserved.
 //
 
 #include "runSolver.hpp"
 #include <assert.h>
 #include "Utility/progress_display.cpp"
 
 void runSolver(Configuration config, Setup &setup){
 
     clock_t tBegin = clock();
     double tMax = config.getTmax();
     // Initialize object
 
     Solution soln = *setup.initialize(config);
     soln.initialize(config);
     MeshBase *mesh = setup.getMesh();
     BcBase *bc = mesh->getBc();
     ModelBase *model = setup.getModel();
     IntBase *integrate = setup.getTimeInt();
     double cfl = config.getCfl();
-    std::vector<std::vector<double> > UNext(mesh->getNcellsTot() , std::vector<double>(mesh->getNelem()));
+    std::vector<std::vector<double> > U(mesh->getNcellsTot() , std::vector<double>(mesh->getNelem()));
     std::vector<double> tInclude = config.getTInclude();
     auto find0 = std::find(tInclude.begin(),tInclude.end(),0);
     if (find0 != tInclude.end()){
         tInclude.erase(find0);
     }
     auto findtMax = std::find(tInclude.begin(),tInclude.end(),tMax);
     if (findtMax == tInclude.end()){
         tInclude.push_back(tMax);
     }
     int nextInclude = 0;
 
     // Initial Data
     double t = 0;
     double tNext;
     double dt;
-    std::vector<std::vector<double>> U = soln.readInitialData();
+    std::vector<std::vector<double>> UNext = soln.readInitialData();
     if (config.getImportSol()){
         soln.importHDF5(t, U);
     }
 
     bc->updateBoundary(U, mesh);
     soln.pushSolution(U, t);
 
     // Run main loops
     int nTInclude = (int) tInclude.size();
 //    setup.exportStencils(U, mesh, t, config);
     while (t < tMax && nextInclude < nTInclude){
-
+        U = UNext;
         dt = cfl*model->calcTimestep(U, mesh);
 //        dt *= 10;
 //        dt = 0.001;
 //        std::cout << dt << "\n";
         if (t+dt >tInclude[nextInclude]){
             dt = tInclude[nextInclude]-t;
             tNext = t + dt;
             nextInclude++;
         }
         else{
             tNext = t + dt;
         }
 
         printProgress(t/tMax);
         integrate->timeInt(UNext, U, t, dt);
 //        if (!isWellDef(UNext)){
 //            clock_t tEnd = clock();
 //            soln.pushSolution(UNext, tNext);
 //            soln.finish((tEnd-tBegin)/CLOCKS_PER_SEC);
 //            std::cout << "sol is not well defined!" << "\n";
 //            assert(isWellDef(UNext));
 //        }
 //        assert(isWellDef(UNext));
         if (model->breaksPositivity(UNext)){
             clock_t tEnd = clock();
             soln.pushSolution(UNext, tNext);
             soln.finish((tEnd-tBegin)/CLOCKS_PER_SEC);
             std::cout << "sol breaks positivity!" << "\n";
             setup.exportStencils(U, mesh, t, config);
             assert(!model->breaksPositivity(UNext));
         }
 //        assert(!model->breaksPositivity(UNext));
         
         soln.pushSolution(UNext, tNext);
-        
-        U = UNext;
+
         t = tNext;
     }
     clock_t tEnd = clock();
+    printResidual(U, UNext);
     soln.finish((tEnd-tBegin)/(double)CLOCKS_PER_SEC);
 };
+
+void printResidual(std::vector<std::vector<double>> &U, std::vector<std::vector<double>> &UNext){
+    int n = U.size();
+    int m = U[0].size();
+    double res = 0;
+    double tmp;
+    for (int i = 0; i < n; i++){
+        for (int j = 0; j < m; j++){
+            tmp = U[i][j]-UNext[i][j];
+            res += tmp*tmp;
+        }
+    }
+
+    std::cout<< "Residual = " << std::sqrt(res/(n*m)) << std::endl;
+}
\ No newline at end of file
diff --git a/ConSol/ConSol/Base/runSolver.hpp b/ConSol/ConSol/Base/runSolver.hpp
index e21079f..bb44e11 100644
--- a/ConSol/ConSol/Base/runSolver.hpp
+++ b/ConSol/ConSol/Base/runSolver.hpp
@@ -1,19 +1,21 @@
 //
 //  runSolver.hpp
 //  ConSol
 //
 //  Created by Fabian Moenkeberg on 30.06.17.
 //  Copyright © 2017 EPFL. All rights reserved.
 //
 
 #ifndef runSolver_hpp
 #define runSolver_hpp
 
 #include <stdio.h>
 //#include "Solution.hpp"
 #include "Setup.hpp"
 #include "../Utility/isWellDef.hpp"
 
 void runSolver(Configuration config, Setup &setup);
 
+void printResidual(std::vector<std::vector<double>> &U, std::vector<std::vector<double>> &UNext);
+
 #endif /* runSolver_hpp */