diff --git a/lyopronto/opt_Pch.py b/lyopronto/opt_Pch.py
index 6a9be783..c272670c 100644
--- a/lyopronto/opt_Pch.py
+++ b/lyopronto/opt_Pch.py
@@ -63,6 +63,10 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
     # Objective function to be minimized to maximize sublimation rate
     def objfun(x):
         return (x[0]-x[4])
+    # Exact gradient of the linear objective, so SLSQP does not
+    # finite-difference it at every point.
+    def objfun_jac(x):
+        return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0])
     # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv]
     x0 = np.array([P0,0.0,Tb0,Tsh0,P0*1.1,Ts0,3.0e-4])    # Initial values
     failures = 0
@@ -71,19 +75,23 @@ def objfun(x):
 
         Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])  # Product resistance [cm^2-hr-Torr/g]
 
-        # Constraints
-        cons = ({'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[0]},  # sublimation front pressure [Torr]
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[1]},    # sublimation rate [kg/hr]
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[2]},    # vial heat transfer balance
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[3]},    # shelf temperature [degC]
-            {'type':'eq','fun':lambda x: x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0])},    # vial heat transfer coefficient [cal/s/K/cm^2]
-            {'type':'eq','fun':lambda x: x[3]-Tsh},    # shelf temperature fixed [degC]
-            {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[0]},  # equipment capability inequlity
-            {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[1]})  # maximum product temperature inequality
+        # Stack the equality constraints into one vector-valued constraint so
+        # SLSQP evaluates and differentiates the whole system once per point
+        # rather than once per component: sublimation front pressure [Torr],
+        # sublimation rate [kg/hr], vial heat transfer balance, shelf
+        # temperature [degC], vial heat transfer coefficient [cal/s/K/cm^2], and fixed shelf temperature [degC]
+        def eq_sys(x, Tsh=Tsh):
+            return np.array(functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)
+                            + (x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0]), x[3]-Tsh))
+        # Inequality constraints: equipment capability and maximum product temperature
+        def ineq_sys(x):
+            return np.array(functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial))
+        cons = ({'type':'eq','fun':eq_sys},
+            {'type':'ineq','fun':ineq_sys})
         # Bounds for the unknowns
         bnds = ((Pchamber['min'],Pchamber.get('max', None)),(0,None),(None,None),(None,None),(0,None),(None,None),(0,None))
         # Minimize the objective function i.e. maximize the sublimation rate
-        res = sp.minimize(objfun,x0,bounds = bnds, constraints = cons)
+        res = sp.minimize(objfun,x0,jac = objfun_jac,bounds = bnds, constraints = cons)
         [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x']    # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2]
         # # Use the results as a guess for the next iteration
         # TODO: decide on appropriate error handling for unsuccessful iterations
diff --git a/lyopronto/opt_Pch_Tsh.py b/lyopronto/opt_Pch_Tsh.py
index 214edf81..23b4ff83 100644
--- a/lyopronto/opt_Pch_Tsh.py
+++ b/lyopronto/opt_Pch_Tsh.py
@@ -52,21 +52,30 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
         Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])  # Product resistance [cm^2-hr-Torr/g]
     
         # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv]
-        def fun(x): 
+        def objfun(x): 
             return x[0]-x[4]    # Objective function to be minimized to maximize sublimation rate
+        # Exact gradient of the linear objective, so SLSQP does not
+        # finite-difference it at every point.
+        def objfun_jac(x):
+            return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0])
         x0 = [P0,0.0,T0,T0,P0,T0,3.0e-4]    # Initial values
-        # Constraints
-        cons = ({'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[0]},  # sublimation front pressure [Torr]
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[1]},    # sublimation rate [kg/hr]
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[2]},    # vial heat transfer balance
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[3]},    # shelf temperature [degC]
-            {'type':'eq','fun':lambda x: x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0])},    # vial heat transfer coefficient [cal/s/K/cm^2]
-            {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[0]},  # equipment capability inequlity
-            {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[1]})  # maximum product temperature inequality
+        # Stack the equality constraints into one vector-valued constraint so
+        # SLSQP evaluates and differentiates the whole system once per point
+        # rather than once per component: sublimation front pressure [Torr],
+        # sublimation rate [kg/hr], vial heat transfer balance, shelf
+        # temperature [degC], vial heat transfer coefficient [cal/s/K/cm^2]
+        def eq_sys(x):
+            return np.array(functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)
+                            + (x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0]),))
+        # Inequality constraints: equipment capability and maximum product temperature
+        def ineq_sys(x):
+            return np.array(functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial))
+        cons = ({'type':'eq','fun':eq_sys},
+            {'type':'ineq','fun':ineq_sys})
         # Bounds for the unknowns
         bnds = ((Pchamber['min'],Pchamber.get('max', None)),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None))
         # Minimize the objective function i.e. maximize the sublimation rate
-        res = sp.minimize(fun,x0,bounds = bnds, constraints = cons)
+        res = sp.minimize(objfun,x0,jac = objfun_jac,bounds = bnds, constraints = cons)
         [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x']    # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2]
 
         # Sublimated ice length
diff --git a/lyopronto/opt_Tsh.py b/lyopronto/opt_Tsh.py
index 20f6b3e9..9a88f71c 100644
--- a/lyopronto/opt_Tsh.py
+++ b/lyopronto/opt_Tsh.py
@@ -59,22 +59,30 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
         Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])  # Product resistance [cm^2-hr-Torr/g]
     
         # Quantities solved for: x = [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv]
-        def fun(x): 
+        def objfun(x): 
             return (x[0]-x[4])    # Objective function to be minimized to maximize sublimation rate
+        # Exact gradient of the linear objective, so SLSQP does not
+        # finite-difference it at every point.
+        def objfun_jac(x):
+            return np.array([1.0,0.0,0.0,0.0,-1.0,0.0,0.0])
         x0 = [Pch,0.0,T0,T0,Pch,T0,3.0e-4]    # Initial values
-        # Constraints
-        cons = ({'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[0]},  # sublimation front pressure [Torr]
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[1]},    # sublimation rate [kg/hr]
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[2]},    # vial heat transfer balance
-            {'type':'eq','fun':lambda x: functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)[3]},    # shelf temperature [degC]
-            {'type':'eq','fun':lambda x: x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0])},    # vial heat transfer coefficient [cal/s/K/cm^2]
-            {'type':'eq','fun':lambda x: x[0]-Pch},    # chamber pressure fixed [Torr]
-            {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[0]},  # equipment capability inequlity
-            {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[1]})  # maximum product temperature inequality
+        # Stack the equality constraints into one vector-valued constraint so
+        # SLSQP evaluates and differentiates the whole system once per point
+        # rather than once per component: sublimation front pressure [Torr],
+        # sublimation rate [kg/hr], vial heat transfer balance, shelf
+        # temperature [degC], vial heat transfer coefficient [cal/s/K/cm^2], and fixed chamber pressure [Torr]
+        def eq_sys(x, Pch=Pch):
+            return np.array(functions.Eq_Constraints(x[0],x[1],x[2],x[3],x[4],x[5],x[6],Lpr0,Lck,vial['Av'],vial['Ap'],Rp)
+                            + (x[6]-functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],x[0]), x[0]-Pch))
+        # Inequality constraints: equipment capability and maximum product temperature
+        def ineq_sys(x):
+            return np.array(functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial))
+        cons = ({'type':'eq','fun':eq_sys},
+            {'type':'ineq','fun':ineq_sys})
         # Bounds for the unknowns
         bnds = ((None,None),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None))
         # Minimize the objective function i.e. maximize the sublimation rate
-        res = sp.minimize(fun,x0,bounds = bnds, constraints = cons)
+        res = sp.minimize(objfun,x0,jac = objfun_jac,bounds = bnds, constraints = cons)
         [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x']    # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2]
 
         # Sublimated ice length