diff --git a/include/cantera/thermo/MixtureFugacityTP.h b/include/cantera/thermo/MixtureFugacityTP.h
index e2b9fc1ed..de9aa1d87 100644
--- a/include/cantera/thermo/MixtureFugacityTP.h
+++ b/include/cantera/thermo/MixtureFugacityTP.h
@@ -551,14 +551,6 @@ protected:
 
     //@}
 
-    class spinodalFunc : public ResidEval
-    {
-    public:
-        spinodalFunc(MixtureFugacityTP* tp);
-        virtual int evalSS(const doublereal t, const doublereal* const y, doublereal* const r);
-        MixtureFugacityTP* m_tp;
-    };
-
 protected:
     virtual void invalidateCache();
 
@@ -604,8 +596,6 @@ protected:
 
     //! Temporary storage for dimensionless reference state entropies
     mutable vector_fp m_s0_R;
-
-    spinodalFunc* fdpdv_;
 };
 }
 
diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp
index 91b70ce93..96eafe94c 100644
--- a/src/thermo/MixtureFugacityTP.cpp
+++ b/src/thermo/MixtureFugacityTP.cpp
@@ -546,20 +546,6 @@ void MixtureFugacityTP::updateMixingExpressions()
 {
 }
 
-MixtureFugacityTP::spinodalFunc::spinodalFunc(MixtureFugacityTP* tp) :
-    m_tp(tp)
-{
-}
-
-int MixtureFugacityTP::spinodalFunc::evalSS(const doublereal t, const doublereal* const y,
-        doublereal* const r)
-{
-    doublereal molarVol = y[0];
-    doublereal pp;
-    r[0] = m_tp->dpdVCalc(m_tp->temperature(), molarVol, pp);
-    return 0;
-}
-
 int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& densLiqGuess,
                              doublereal& densGasGuess, doublereal& liqGRT, doublereal& gasGRT)
 {
diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp
index 67d850778..73a739542 100644
--- a/src/thermo/RedlichKwongMFTP.cpp
+++ b/src/thermo/RedlichKwongMFTP.cpp
@@ -6,11 +6,13 @@
 #include "cantera/thermo/RedlichKwongMFTP.h"
 #include "cantera/thermo/mix_defs.h"
 #include "cantera/thermo/ThermoFactory.h"
-#include "cantera/numerics/RootFind.h"
 #include "cantera/base/stringUtils.h"
 #include "cantera/base/ctml.h"
 
+#include <boost/math/tools/roots.hpp>
+
 using namespace std;
+namespace bmt = boost::math::tools;
 
 namespace Cantera
 {
@@ -899,46 +901,46 @@ doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa,
 
 doublereal RedlichKwongMFTP::densSpinodalLiquid() const
 {
-    if (NSolns_ != 3) {
+    double Vroot[3];
+    double T = temperature();
+    int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
+    if (nsol != 3) {
         return critDensity();
     }
-    double vmax = Vroot_[1];
-    double vmin = Vroot_[0];
-    RootFind rf(fdpdv_);
-    rf.setPrintLvl(10);
-    rf.setTol(1.0E-5, 1.0E-10);
-    rf.setFuncIsGenerallyDecreasing(true);
 
-    double vbest = 0.5 * (Vroot_[0]+Vroot_[1]);
-    double funcNeeded = 0.0;
-    int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
-    if (status != ROOTFIND_SUCCESS) {
-        throw CanteraError("  RedlichKwongMFTP::densSpinodalLiquid() ", "didn't converge");
-    }
+    auto resid = [this, T](double v) {
+        double pp;
+        return dpdVCalc(T, v, pp);
+    };
+
+    boost::uintmax_t maxiter = 100;
+    std::pair<double, double> vv = bmt::toms748_solve(
+        resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
+
     doublereal mmw = meanMolecularWeight();
-    return mmw / vbest;
+    return mmw / (0.5 * (vv.first + vv.second));
 }
 
 doublereal RedlichKwongMFTP::densSpinodalGas() const
 {
-    if (NSolns_ != 3) {
+    double Vroot[3];
+    double T = temperature();
+    int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
+    if (nsol != 3) {
         return critDensity();
     }
-    double vmax = Vroot_[2];
-    double vmin = Vroot_[1];
-    RootFind rf(fdpdv_);
-    rf.setPrintLvl(10);
-    rf.setTol(1.0E-5, 1.0E-10);
-    rf.setFuncIsGenerallyIncreasing(true);
 
-    double vbest = 0.5 * (Vroot_[1]+Vroot_[2]);
-    double funcNeeded = 0.0;
-    int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
-    if (status != ROOTFIND_SUCCESS) {
-        throw CanteraError(" RedlichKwongMFTP::densSpinodalGas() ", "didn't converge");
-    }
+    auto resid = [this, T](double v) {
+        double pp;
+        return dpdVCalc(T, v, pp);
+    };
+
+    boost::uintmax_t maxiter = 100;
+    std::pair<double, double> vv = bmt::toms748_solve(
+        resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
+
     doublereal mmw = meanMolecularWeight();
-    return mmw / vbest;
+    return mmw / (0.5 * (vv.first + vv.second));
 }
 
 doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const