diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index a73b13c73e310c805405a2319e4260a236e083cf..45574e241f36afcbda47f9c92d5c186cefe2fc7e 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -323,6 +323,33 @@ int main(int argc, char *argv[]) {
                         << std::endl;
         }
 
+        // Comparison with the analytic solution of a velocity stepping test
+        // with the Ruina state evolution law.
+        // Jumps at 120 and 360 timesteps; v1 = .6 * v2;
+        if (parset.get<bool>("printVelocitySteppingComparison")) {
+          double const v = u4_diff[first_frictional_node].two_norm() / h / L;
+
+          if (run >= 120) {
+            double const euler = (*s4_new)[first_frictional_node];
+            double direct;
+            if (run < 360) {
+              double const v2 = v;
+              double const v1 = 0.6 * v2;
+              direct = std::log(
+                  1.0 / v2 *
+                  std::pow((v2 / v1), std::exp(-v2 * (run - 120) * h)));
+            } else {
+              double const v1 = v;
+              double const v2 = v1 / 0.6;
+              direct = std::log(
+                  1.0 / v1 *
+                  std::pow((v1 / v2), std::exp(-v1 * (run - 360) * h)));
+            }
+            std::cout << "[" << run << "]: " << euler << " " << direct << " "
+                      << (euler - direct) / direct << std::endl;
+          }
+        }
+
         // Record the coefficient of friction at a fixed node
         if (parset.get<bool>("printCoefficient")) {
           double const V = u4_diff[first_frictional_node].two_norm();
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index 62fb09c45cf03a4c2c94aac126094dd282718efe..58fa7b7968e0fbe8053b492eca340e172e5447c7 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -9,6 +9,8 @@ printProgress = true
 writeEvolution = false
 writeVTK = false
 
+printVelocitySteppingComparison = false
+
 [grid]
 refinements = 4