diff --git a/include/cantera/tpx/Sub.h b/include/cantera/tpx/Sub.h
index 6a4a96512..7c019c945 100644
--- a/include/cantera/tpx/Sub.h
+++ b/include/cantera/tpx/Sub.h
@@ -128,35 +128,41 @@ public:
     //! Specific heat at constant volume [J/kg/K]
     virtual double cv() {
         double Tsave = T, dt = 1.e-4*T;
-        set_T(Tsave - dt);
+        double T1 = std::max(Tmin(), Tsave - dt);
+        double T2 = std::min(Tmax(), Tsave + dt);
+        set_T(T1);
         double s1 = s();
-        set_T(Tsave + dt);
+        set_T(T2);
         double s2 = s();
         set_T(Tsave);
-        return T*(s2 - s1)/(2.0*dt);
+        return T*(s2 - s1)/(T2-T1);
     }
 
     //! Specific heat at constant pressure [J/kg/K]
     virtual double cp() {
         double Tsave = T, dt = 1.e-4*T;
+        double T1 = std::max(Tmin(), Tsave - dt);
+        double T2 = std::min(Tmax(), Tsave + dt);
         double p0 = P();
-        Set(PropertyPair::TP, Tsave - dt, p0);
+        Set(PropertyPair::TP, T1, p0);
         double s1 = s();
-        Set(PropertyPair::TP, Tsave + dt, p0);
+        Set(PropertyPair::TP, T2, p0);
         double s2 = s();
         Set(PropertyPair::TP, Tsave, p0);
-        return T*(s2 - s1)/(2.0*dt);
+        return T*(s2 - s1)/(T2-T1);
     }
 
     virtual double thermalExpansionCoeff() {
         double Tsave = T, dt = 1.e-4*T;
+        double T1 = std::max(Tmin(), Tsave - dt);
+        double T2 = std::min(Tmax(), Tsave + dt);
         double p0 = P();
-        Set(PropertyPair::TP, Tsave - dt, p0);
+        Set(PropertyPair::TP, T1, p0);
         double v1 = v();
-        Set(PropertyPair::TP, Tsave + dt, p0);
+        Set(PropertyPair::TP, T2, p0);
         double v2 = v();
         Set(PropertyPair::TP, Tsave, p0);
-        return (v2 - v1)/((v2 + v1)*dt);
+        return (v2 - v1)/((v2 + v1)*(T2-T1));
     }
 
     virtual double isothermalCompressibility() {
diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index c54210a34..5ddc6e842 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -18,3 +18,37 @@ class TestPureFluid(utilities.CanteraTest):
         self.water.TX = 500, 0.8
         self.assertNear(self.water.T, 500)
         self.assertNear(self.water.X, 0.8)
+
+    def test_set_minmax(self):
+        self.water.TP = self.water.min_temp, 101325
+        self.assertNear(self.water.T, self.water.min_temp)
+
+        self.water.TP = self.water.max_temp, 101325
+        self.assertNear(self.water.T, self.water.max_temp)
+
+    def check_fd_properties(self, T1, P1, T2, P2, tol):
+        # Properties which are computed as finite differences
+        self.water.TP = T1, P1
+        cp1 = self.water.cp_mass
+        cv1 = self.water.cv_mass
+        k1 = self.water.isothermal_compressibility
+        alpha1 = self.water.thermal_expansion_coeff
+
+        self.water.TP = T2, P2
+        cp2 = self.water.cp_mass
+        cv2 = self.water.cv_mass
+        k2 = self.water.isothermal_compressibility
+        alpha2 = self.water.thermal_expansion_coeff
+
+        self.assertNear(cp1, cp2, tol)
+        self.assertNear(cv1, cv2, tol)
+        self.assertNear(k1, k2, tol)
+        self.assertNear(alpha1, alpha2, tol)
+
+    def test_properties_near_min(self):
+        self.check_fd_properties(self.water.min_temp*(1+1e-5), 101325,
+                              self.water.min_temp*(1+1e-4), 101325, 1e-2)
+
+    def test_properties_near_max(self):
+        self.check_fd_properties(self.water.max_temp*(1-1e-5), 101325,
+                              self.water.max_temp*(1-1e-4), 101325, 1e-2)
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 9fbbfa0f1..9167ed4ef 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -5,6 +5,7 @@
  */
 #include "cantera/tpx/Sub.h"
 #include "cantera/base/stringUtils.h"
+#include "cantera/base/global.h"
 
 using std::string;
 using namespace Cantera;
@@ -35,9 +36,9 @@ double Substance::dPsdT()
 {
     double tsave = T;
     double ps1 = Ps();
-    set_T(T + DeltaT);
+    T = T + DeltaT;
     double dpdt = (Ps() - ps1)/DeltaT;
-    set_T(tsave);
+    T = tsave;
     return dpdt;
 }
 
@@ -387,7 +388,7 @@ int Substance::Lever(int itp, double sat, double val, propertyFlag::type ifunc)
         if (sat >= Tcrit()) {
             return 0;
         }
-        set_T(sat);
+        T = sat;
         psat = Ps();
     } else if (itp == Pgiven) {
         if (sat >= Pcrit()) {
@@ -506,11 +507,7 @@ void Substance::set_xy(propertyFlag::type ifx, propertyFlag::type ify,
         }
         v_here += dv;
         t_here += dt;
-        if (t_here >= Tmax()) {
-            t_here = Tmax() - 0.001;
-        } else if (t_here <= Tmin()) {
-            t_here = Tmin() + 0.001;
-        }
+        t_here = clip(t_here, Tmin(), Tmax());
         if (v_here <= 0.0) {
             v_here = 0.0001;
         }
@@ -580,7 +577,7 @@ void Substance::set_TPp(double Temp, double Pressure)
     int LoopCount = 0;
 
     double v_save = 1.0/Rho;
-    set_T(Temp);
+    T = Temp;
     v_here = vp();
 
     // loop