diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index 6b9862331448b9d3bcb7c4ae8a4e1f21179cf513..5432a54e22a1400efbb31118c6328b85fd31d980 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -299,7 +299,7 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
         return defaultdict(lambda: [],
                 ((irow, [
                     (icol,
-                        coeff.subs(self._rscale_symbol, rscale)
+                        coeff.xreplace({self._rscale_symbol: rscale})
                         if isinstance(coeff, sp.Basic)
                         else coeff)
                     for icol, coeff in row])
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 0ddbb428108ff0bd52dd896f88ba1115a66428c5..4e226cd96edfd7c72973116913a11474775b694c 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -25,7 +25,7 @@ THE SOFTWARE.
 from six.moves import range, zip
 import sumpy.symbolic as sym  # noqa
 
-from sumpy.symbolic import vector_subs
+from sumpy.symbolic import vector_xreplace
 from sumpy.expansion import (
     ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion,
     HelmholtzConformingVolumeTaylorExpansion)
@@ -105,7 +105,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
         if self.kernel.has_efficient_scale_adjustment:
             return (
                     self.kernel.adjust_for_kernel_scaling(
-                        vector_subs(
+                        vector_xreplace(
                             expr,
                             bvec, bvec * rscale**-1),
                         rscale, nderivatives)
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index ddfba6829ca6c1fdb4bfc6930cc9352e7c1d89e0..15ad22a522981d196e6dcaccb81ed22a4f0d2d7f 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -340,7 +340,7 @@ class ExpressionKernel(Kernel):
             raise ValueError("dist_vec length does not match expected dimension")
 
         from sumpy.symbolic import Symbol
-        expr = expr.subs(dict(
+        expr = expr.xreplace(dict(
             (Symbol("d%d" % i), dist_vec_i)
             for i, dist_vec_i in enumerate(scaled_dist_vec)
             ))
diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py
index 49fa405929ed713579b27e2b52d6af2493da482b..c4cf15c9d0f6232943e330b16932fd95eed1f96a 100644
--- a/sumpy/symbolic.py
+++ b/sumpy/symbolic.py
@@ -191,14 +191,14 @@ def make_sym_vector(name, components):
             [sym.Symbol("%s%d" % (name, i)) for i in range(components)])
 
 
-def vector_subs(expr, from_vec, to_vec):
-    subs_pairs = []
+def vector_xreplace(expr, from_vec, to_vec):
+    substs = {}
     assert (from_vec.rows, from_vec.cols) == (to_vec.rows, to_vec.cols)
     for irow in range(from_vec.rows):
         for icol in range(from_vec.cols):
-            subs_pairs.append((from_vec[irow, icol], to_vec[irow, icol]))
+            substs[from_vec[irow, icol]] = to_vec[irow, icol]
 
-    return expr.subs(dict(subs_pairs))
+    return expr.xreplace(substs)
 
 
 def find_power_of(base, prod):