Handle edge case in AdaptiveIntegrator::step().

1 files changed, 10 insertions, 6 deletions

diff --git a/src/vZ/Adaptive.hpp b/src/vZ/Adaptive.hpp
index 256bda6..05fc439 100644
--- a/src/vZ/Adaptive.hpp
+++ b/src/vZ/Adaptive.hpp
@@ -76,30 +76,34 @@ namespace vZ
   GenericAdaptiveIntegrator<Y>::step()
   {
     static const Scalar S = Scalar(19)/Scalar(20); // Arbitrary saftey factor
-    Scalar newH;
+    Scalar newH = this->h();
     Y y;
 
     // Attempt the integration step in a loop
-    bool rejected = true;
-    while (rejected) {
+    while (true) {
       KVector k = calculateK(m_a);
       y         = calculateY(k, m_b);
       Y yStar   = calculateY(k, m_bStar);
 
       // Get an error estimate
+
       using std::abs;
       using std::pow;
       Scalar delta = abs(y - yStar);
-      Scalar scale = m_atol + std::max(abs(y), abs(this->y()))*m_rtol;
 
-      newH = S*this->h()*pow(scale/delta, Scalar(1)/m_order);
+      if (delta == Scalar(0)) {
+        break;
+      }
+
+      Scalar scale = m_atol + std::max(abs(y), abs(this->y()))*m_rtol;
+      newH         = S*this->h()*pow(scale/delta, Scalar(1)/m_order);
 
       if (delta > scale) {
         // Reject the step
         this->h(newH);
         ++m_rejections;
       } else {
-        rejected = false;
+        break;
       }
     }