diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp
index 8bf88ed7c..0ab925641 100644
--- a/src/oneD/MultiNewton.cpp
+++ b/src/oneD/MultiNewton.cpp
@@ -26,17 +26,133 @@ using namespace std;
 namespace Cantera
 {
 
-//----------------------------------------------------------
-//                  function declarations
-//----------------------------------------------------------
+// unnamed-namespace for local helpers
+namespace {
 
-// declarations for functions in newton_utils.h
-doublereal bound_step(const doublereal* x,
-                      const doublereal* step, Domain1D& r, int loglevel=0);
+class Indx
+{
+public:
+    Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
+    size_t m_nv, m_np;
+    size_t operator()(size_t m, size_t j) {
+        return j*m_nv + m;
+    }
+};
+
+
+/**
+ * Return a damping coefficient that keeps the solution after taking one
+ * Newton step between specified lower and upper bounds. This function only
+ * considers one domain.
+ */
+doublereal bound_step(const doublereal* x, const doublereal* step,
+                      Domain1D& r, int loglevel)
+{
+
+    char buf[100];
+    size_t np = r.nPoints();
+    size_t nv = r.nComponents();
+    Indx index(nv, np);
+    doublereal above, below, val, newval;
+    size_t m, j;
+    doublereal fbound = 1.0;
+    bool wroteTitle = false;
+    for (m = 0; m < nv; m++) {
+        above = r.upperBound(m);
+        below = r.lowerBound(m);
+
+        for (j = 0; j < np; j++) {
+            val = x[index(m,j)];
+            if (loglevel > 0) {
+                if (val > above + 1.0e-12 || val < below - 1.0e-12) {
+                    sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n",
+                            int2str(r.domainIndex()).c_str(),
+                            r.componentName(m).c_str(), int2str(j).c_str(),
+                            val, below, above);
+                    writelog(string("\nERROR: solution out of bounds.\n")+buf);
+                }
+            }
+
+            newval = val + step[index(m,j)];
+
+            if (newval > above) {
+                fbound = std::max(0.0, std::min(fbound,
+                                                (above - val)/(newval - val)));
+            } else if (newval < below) {
+                fbound = std::min(fbound, (val - below)/(val - newval));
+            }
+
+            if (loglevel > 1 && (newval > above || newval < below)) {
+                if (!wroteTitle) {
+                    writelog("\nNewton step takes solution out of bounds.\n\n");
+                    sprintf(buf,"  %12s  %12s  %4s  %10s  %10s  %10s  %10s\n",
+                            "domain","component","pt","value","step","min","max");
+                    wroteTitle = true;
+                    writelog(buf);
+                }
+                sprintf(buf, "          %4s  %12s  %4s  %10.3e  %10.3e  %10.3e  %10.3e\n",
+                        int2str(r.domainIndex()).c_str(),
+                        r.componentName(m).c_str(), int2str(j).c_str(),
+                        val, step[index(m,j)], below, above);
+                writelog(buf);
+            }
+        }
+    }
+    return fbound;
+}
+
+
+/**
+ * This function computes the square of a weighted norm of a step
+ * vector for one domain.
+ *
+ * @param x     Solution vector for this domain.
+ * @param step  Newton step vector for this domain.
+ * @param r     Object representing the domain. Used to get tolerances,
+ *              number of components, and number of points.
+ *
+ * The return value is
+ * \f[
+ *    \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
+ * \f]
+ * where the error weight for solution component \f$n\f$ is given by
+ * \f[
+ *     w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
+ * \f]
+ * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for
+ * component n, and multiplies the average magnitude of
+ * solution component n in the domain. The second term,
+ * \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
+ * n.
+ *
+ */
 doublereal norm_square(const doublereal* x,
-                       const doublereal* step, Domain1D& r);
+                       const doublereal* step, Domain1D& r)
+{
+    doublereal f, ewt, esum, sum = 0.0;
+    size_t n, j;
+    doublereal f2max = 0.0;
+    size_t nv = r.nComponents();
+    size_t np = r.nPoints();
 
+    for (n = 0; n < nv; n++) {
+        esum = 0.0;
+        for (j = 0; j < np; j++) {
+            esum += fabs(x[nv*j + n]);
+        }
+        ewt = r.rtol(n)*esum/np + r.atol(n);
+        for (j = 0; j < np; j++) {
+            f = step[nv*j + n]/ewt;
+            sum += f*f;
+            if (f*f > f2max) {
+                f2max = f*f;
+            }
+        }
+    }
+    return sum;
+}
 
+} // end unnamed-namespace
 
 //-----------------------------------------------------------
 //                  constants
@@ -424,7 +540,5 @@ void MultiNewton::releaseWorkArray(doublereal* work)
 {
     m_workarrays.push_back(work);
 }
-}
 
-
-// $Log: Newton.cpp,v
+} // end namespace Cantera
diff --git a/src/oneD/newton_utils.cpp b/src/oneD/newton_utils.cpp
deleted file mode 100644
index a3ef61812..000000000
--- a/src/oneD/newton_utils.cpp
+++ /dev/null
@@ -1,138 +0,0 @@
-/**
- *  @file newton_utils.cpp
- */
-
-#include "cantera/base/ct_defs.h"
-#include "cantera/oneD/Domain1D.h"
-
-#include <cstdio>
-
-using namespace std;
-
-namespace Cantera
-{
-
-class Indx
-{
-public:
-    Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
-    size_t m_nv, m_np;
-    size_t operator()(size_t m, size_t j) {
-        return j*m_nv + m;
-    }
-};
-
-
-/**
- * Return a damping coefficient that keeps the solution after taking one
- * Newton step between specified lower and upper bounds. This function only
- * considers one domain.
- */
-doublereal bound_step(const doublereal* x, const doublereal* step,
-                      Domain1D& r, int loglevel)
-{
-
-    char buf[100];
-    size_t np = r.nPoints();
-    size_t nv = r.nComponents();
-    Indx index(nv, np);
-    doublereal above, below, val, newval;
-    size_t m, j;
-    doublereal fbound = 1.0;
-    bool wroteTitle = false;
-    for (m = 0; m < nv; m++) {
-        above = r.upperBound(m);
-        below = r.lowerBound(m);
-
-        for (j = 0; j < np; j++) {
-            val = x[index(m,j)];
-            if (loglevel > 0) {
-                if (val > above + 1.0e-12 || val < below - 1.0e-12) {
-                    sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n",
-                            int2str(r.domainIndex()).c_str(),
-                            r.componentName(m).c_str(), int2str(j).c_str(),
-                            val, below, above);
-                    writelog(string("\nERROR: solution out of bounds.\n")+buf);
-                }
-            }
-
-            newval = val + step[index(m,j)];
-
-            if (newval > above) {
-                fbound = std::max(0.0, std::min(fbound,
-                                                (above - val)/(newval - val)));
-            } else if (newval < below) {
-                fbound = std::min(fbound, (val - below)/(val - newval));
-            }
-
-            if (loglevel > 1 && (newval > above || newval < below)) {
-                if (!wroteTitle) {
-                    writelog("\nNewton step takes solution out of bounds.\n\n");
-                    sprintf(buf,"  %12s  %12s  %4s  %10s  %10s  %10s  %10s\n",
-                            "domain","component","pt","value","step","min","max");
-                    wroteTitle = true;
-                    writelog(buf);
-                }
-                sprintf(buf, "          %4s  %12s  %4s  %10.3e  %10.3e  %10.3e  %10.3e\n",
-                        int2str(r.domainIndex()).c_str(),
-                        r.componentName(m).c_str(), int2str(j).c_str(),
-                        val, step[index(m,j)], below, above);
-                writelog(buf);
-            }
-        }
-    }
-    return fbound;
-}
-
-
-
-/**
- * This function computes the square of a weighted norm of a step
- * vector for one domain.
- *
- * @param x     Solution vector for this domain.
- * @param step  Newton step vector for this domain.
- * @param r     Object representing the domain. Used to get tolerances,
- *              number of components, and number of points.
- *
- * The return value is
- * \f[
- *    \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
- * \f]
- * where the error weight for solution component \f$n\f$ is given by
- * \f[
- *     w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
- * \f]
- * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for
- * component n, and multiplies the average magnitude of
- * solution component n in the domain. The second term,
- * \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
- * n.
- *
- */
-doublereal norm_square(const doublereal* x,
-                       const doublereal* step, Domain1D& r)
-{
-    doublereal f, ewt, esum, sum = 0.0;
-    size_t n, j;
-    doublereal f2max = 0.0;
-    size_t nv = r.nComponents();
-    size_t np = r.nPoints();
-
-    for (n = 0; n < nv; n++) {
-        esum = 0.0;
-        for (j = 0; j < np; j++) {
-            esum += fabs(x[nv*j + n]);
-        }
-        ewt = r.rtol(n)*esum/np + r.atol(n);
-        for (j = 0; j < np; j++) {
-            f = step[nv*j + n]/ewt;
-            sum += f*f;
-            if (f*f > f2max) {
-                f2max = f*f;
-            }
-        }
-    }
-    return sum;
-}
-}