diff --git a/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java b/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java
index b71439966..d7fc189f8 100644
--- a/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java
+++ b/src/main/java/org/apache/commons/math4/analysis/solvers/MullerSolver.java
@@ -29,8 +29,6 @@ import org.apache.commons.math4.util.FastMath;
  * <p>
  * Muller's method applies to both real and complex functions, but here we
  * restrict ourselves to real functions.
- * This class differs from {@link MullerSolver} in the way it avoids complex
- * operations.</p><p>
  * Muller's original method would have function evaluation at complex point.
  * Since our f(x) is real, we have to find ways to avoid that. Bracketing
  * condition is one way to go: by requiring bracketing in every iteration,
@@ -161,6 +159,13 @@ public class MullerSolver extends AbstractUnivariateSolver {
             // xplus and xminus are two roots of parabola and at least
             // one of them should lie in (x0, x2)
             final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
+
+            // XXX debug
+            if (!isSequence(x0, x, x2)) {
+                System.out.println("x=" + x + " x0=" + x0 + " x2=" + x2);
+                throw new org.apache.commons.math4.exception.MathInternalError();
+            }
+
             final double y = computeObjectiveValue(x);
 
             // check for convergence
diff --git a/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java b/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java
index e82be5230..2bb0b1bdf 100644
--- a/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java
+++ b/src/test/java/org/apache/commons/math4/analysis/solvers/MullerSolverTest.java
@@ -147,4 +147,35 @@ public final class MullerSolverTest {
             // expected
         }
     }
+
+    @Test
+    public void testMath1333() {
+        final UnivariateFunction logFunction = new UnivariateFunction() {
+                private double log1pe(double x) {
+                    if (x > 0) {
+                        return x + FastMath.log1p(FastMath.exp(-x));
+                    } else {
+                        return FastMath.log1p(FastMath.exp(x));
+                    }
+                }
+
+                @Override
+                public double value(double x) {
+                    final double a = 0.15076136473214652;
+                    final double b = 4.880819340168248;
+                    final double c = -2330.4196672490493;
+                    final double d = 1.1871451743330544E-16;
+                    //aa*log(1+e^(bbx+c))+d - 0.01 * x - 20 * 0.01
+                    return a * a * log1pe(b * b * x + c) + d - 0.01 * x - 20 * 0.01;
+                }
+            };
+
+        final UnivariateSolver solver = new MullerSolver(0.25);
+        final double min = 20;
+        final double max = 100.04173804515072;
+        final double result = solver.solve(1000, logFunction, min, max, 100 / (double) 3);
+
+        Assert.assertTrue(result + " < " + min, result >= min);
+        Assert.assertTrue(result + " > " + max, result <= max);
+    }
 }