diff --git a/dev/Tickets/882 - tabular phase spec.py b/dev/Tickets/882 - tabular phase spec.py
index 9fd4a2ad..14ae2980 100644
--- a/dev/Tickets/882 - tabular phase spec.py	
+++ b/dev/Tickets/882 - tabular phase spec.py	
@@ -1,8 +1,11 @@
+# coding: utf-8
+
 import msgpack, zlib, StringIO
 import matplotlib.pyplot as plt, numpy as np, CoolProp
 
 # Modify path as necessary
-with open(r'C:\Users\Belli\.CoolProp\Tables\HelmholtzEOSBackend(Water[1.0000000000])/single_phase_logpT.bin.z','rb') as fp:
+user_paths = [r'C:\Users\Belli',r'C:\Users\jowr']
+with open(user_paths[1]+r'\.CoolProp\Tables\HelmholtzEOSBackend(Water[1.0000000000])/single_phase_logpT.bin.z','rb') as fp:
     ph = zlib.decompress(fp.read())
     values = msgpack.load(StringIO.StringIO(ph))
     revision, matrices = values[0:2]
@@ -12,9 +15,11 @@ with open(r'C:\Users\Belli\.CoolProp\Tables\HelmholtzEOSBackend(Water[1.00000000
     Ts = np.linspace(Tt, Tc)
     ps = CoolProp.CoolProp.PropsSI('P','T',Ts,'Q',0,'Water')
     
-    plt.plot(T, p, '.')
+    plt.plot(T, p, '.', color='gray')
     plt.plot(Ts, ps, 'k', lw = 2)
     plt.yscale('log')
+    plt.xlabel('Temperature / K')
+    plt.ylabel('Pressure / Pa')
     plt.show()
     
     
@@ -28,30 +33,45 @@ print AS.rhomolar_critical()
 pt = AS.keyed_output(CoolProp.iP_triple)
 pc = AS.p_critical()
 
+verbose = False
+dTs = np.power(10.0,-np.arange(-1,4))
+
 for p in np.logspace(np.log10(pt*1.05),np.log10(pc*0.95)):
     AS.update(CoolProp.PQ_INPUTS, p, 0)
     rhoL = AS.rhomolar()
     Ts = AS.T()    
     # Liquid side
-    for specify_phase in [True]:
+    for specify_phase in [False,True]:
         if specify_phase:
             AS.specify_phase(CoolProp.iphase_liquid)
         else:
             AS.unspecify_phase()
-        for dT in [10,1,0.01,0.001]:
-            print p, Ts-dT
+        for dT in dTs:
+            if verbose: print p, Ts-dT
             try:
                 AS.update(CoolProp.PT_INPUTS, p, Ts-dT)
-                print p, Ts-dT, AS.rhomolar()
+                if verbose: print p, Ts-dT, AS.rhomolar()
             except BaseException as BE:
-                print 'ERROR:', BE
-        print p, Ts, rhoL
-    print ' '
+                if specify_phase: print 'Liquid error:', BE
+                else:             print 'Liquid issue:', BE
+        if verbose: print p, Ts, rhoL
+    # Gaseous side
+    for specify_phase in [False,True]:
+        if specify_phase:
+            AS.specify_phase(CoolProp.iphase_gas)
+        else:
+            AS.unspecify_phase()
+        for dT in dTs:
+            if verbose: print p, Ts+dT
+            try:
+                AS.update(CoolProp.PT_INPUTS, p, Ts+dT)
+                if verbose: print p, Ts+dT, AS.rhomolar()
+            except BaseException as BE:
+                if specify_phase: print '   Gas error:', BE
+                else:             print '   Gas issue:', BE
+        if verbose: print p, Ts, rhoL
     
-                
-
-
-# In[ ]:
-
-
-
+    
+    
+    
+    
\ No newline at end of file
diff --git a/dev/Tickets/883 - CO2 and water stability.py b/dev/Tickets/883 - CO2 and water stability.py
new file mode 100644
index 00000000..c4a7a695
--- /dev/null
+++ b/dev/Tickets/883 - CO2 and water stability.py	
@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from __future__ import print_function, division
+import numpy as np
+import CoolProp 
+from CoolProp.CoolProp import PropsSI, AbstractState
+import sys
+
+back = "HEOS"
+fluids = ["water", "CO2"]
+
+num_p = 10
+num_h = 500
+
+msgs = ""
+
+for fluid in fluids:
+    state = AbstractState(back, fluid)    
+    fluid = back+"::"+fluid
+    p_lo = 1.05*PropsSI("pmin","",0,"",0,fluid)
+    p_hi = 0.95*PropsSI("pcrit","",0,"",0,fluid)
+    p_r = np.logspace(np.log10(p_lo), np.log10(p_hi), num_p)
+    h_r = []; d_r = []; d_o = []
+    for p in p_r:
+        if state.has_melting_line():
+            T_m = state.melting_line(CoolProp.constants.iT, CoolProp.constants.iP, p)
+            h_m = -np.Inf
+            dT = 0.0
+            while dT < 2.0 and not np.isfinite(h_m):
+                dT += 0.1
+                T_t = T_m + dT 
+                try: 
+                    h_m = PropsSI("H", "P", p, "T", T_t, fluid)
+                except:
+                    pass
+        else:
+            h_m = -np.Inf
+        h1 = PropsSI("H", "P", p, "Q", 0, fluid)
+        dh = 0.15*(PropsSI("H", "P", p, "Q", 1, fluid)-h1)
+        h0 = np.maximum(h1 - dh,h_m) 
+        h = np.linspace(h0,h1,num=num_h)
+        err = [""]*h.size
+        d = np.empty_like(h)
+        for i,h_i in enumerate(h):
+            try: 
+                d[i] = PropsSI("D", "P", p, "H", h_i, fluid)
+                err[i] = ""
+            except Exception as e:
+                d[i] = np.NaN
+                err[i] = str(e)
+
+        err = np.array(err)
+        valid = np.isfinite(d)
+        invalid = np.logical_not(valid)
+        for i in range(np.sum(invalid)):
+            msgs += "[Error] {0} at h={1} J/kg, p={2} Pa: {3}\n".format(fluid,h[invalid][i],p,err[invalid][i])
+
+        print("[Error] {0} at p={1:8.4f} bar: Failed for {2:4.1f}% of the calls.".format(fluid, p/1e5, np.sum(invalid)/h.size*1e2))
+
+#print(msgs)
+
+sys.exit()
\ No newline at end of file