diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index cb89749a4a08fa3ff5d2fdb24e8344c8d9bc1e65..0e40b35f2be3eaa6de583985213b8dc2a7462777 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -138,10 +138,8 @@ int main (int argc, char *argv[]) try
 
     rodX.back().q = Quaternion<double>(axis, M_PI*angle/180);
 
-    std::cout << "Right boundary orientation:" << std::endl;
-    std::cout << "director 0:  " << rodX[rodX.size()-1].q.director(0) << std::endl;
-    std::cout << "director 1:  " << rodX[rodX.size()-1].q.director(1) << std::endl;
-    std::cout << "director 2:  " << rodX[rodX.size()-1].q.director(2) << std::endl;
+    // Backup initial rod iterate for later reference
+    RodSolutionType initialIterateRod = rodX;
 
     int toplevel = rodGrid.maxLevel();
 
@@ -467,7 +465,9 @@ int main (int argc, char *argv[]) try
     // from zero anyways
     oldError += computeEnergyNormSquared(exactSol3d, *hessian3d);
     
-    /** \todo Rod error still missing */
+    // Error of the initial rod iterate
+    RodDifferenceType rodDifference = computeRodDifference(initialIterateRod, exactSolRod);
+    oldError += computeEnergyNormSquared(rodDifference, hessianRod);
 
     oldError = std::sqrt(oldError);