diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 6fbc95c173549a5a73470dae08e0a95b5c8e0583..382c82b1996db1c8a50b41abbb8fc057a55b8717 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -315,7 +315,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
     solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u,
                        parset.sub("u0.solver"));
 
-    //print(u, "initial u:");
+    print(u, "initial u:");
 
     // Initial acceleration: Computed in agreement with Ma = ell0 - Au
     // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
@@ -348,7 +348,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
     solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a,
                        parset.sub("a0.solver"));
 
-    //print(a, "initial a:");
+    print(a, "initial a:");
   }
 
 private:
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index 4cd6024a74fd2904dab51b878371e24b1a01e3bb..c14b16f154a0ca52d5474cbc1e0e19fc2feee57c 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -118,7 +118,8 @@
    void print(const std::vector<T>& x, std::string message){
       std::cout << message << std::endl;
       for (size_t i=0; i<x.size(); i++) {
-          std::cout << x[i] << std::endl;
+          //std::cout << x[i] << std::endl;
+          print(x[i], "");
       }
       std::cout << std::endl << std::endl;
    }