From c3f2867872cbabf04a52381449327527180762dd Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.h.bell@gmail.com>
Date: Wed, 6 Apr 2016 22:01:56 -0600
Subject: [PATCH] Secant behaves properly now when y1 approx. y2; closes #1015

---
 src/Solvers.cpp | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/src/Solvers.cpp b/src/Solvers.cpp
index db33770804..cd77791339 100644
--- a/src/Solvers.cpp
+++ b/src/Solvers.cpp
@@ -152,13 +152,14 @@ double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter
     x = x0;
     while (iter < 2 || std::abs(fval) > ftol)
     {
+        if (f->input_not_in_range(x)){
+            throw ValueError(format("Input [%g] is out of range",x));
+        }
+        
         fval = f->call(x);
         dfdx = f->deriv(x);
         d2fdx2 = f->second_deriv(x);
         
-        if (f->input_not_in_range(x)){
-            throw ValueError(format("Input [%g] is out of range",x));
-        }
         if (!ValidNumber(fval)){
             throw ValueError("Residual function in Halley returned invalid number");
         };
@@ -168,8 +169,6 @@ double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter
         
         dx = -(2*fval*dfdx)/(2*POW2(dfdx)-fval*d2fdx2);
 
-        
-
         x += dx;
 
         if (std::abs(dx/x) < 10*DBL_EPSILON){
@@ -212,6 +211,10 @@ double Secant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter)
         if (iter==1){x1=x0; x=x1;}
         if (iter==2){x2=x0+dx; x=x2;}
         if (iter>2) {x=x2;}
+        
+            if (f->input_not_in_range(x)){
+                throw ValueError(format("Input [%g] is out of range",x));
+            }
 
             fval = f->call(x);
 
@@ -231,6 +234,10 @@ double Secant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter)
                 return x;
             }
             y2=fval;
+            double deltay = y2-y1;
+            if (iter > 2 && std::abs(deltay)<1e-14){
+                return x;
+            }
             x3=x2-y2/(y2-y1)*(x2-x1);
             y1=y2; x1=x2; x2=x3;