diff --git a/benchmarks/bench_translations.py b/benchmarks/bench_translations.py
index 8d6cfdd882be3c8bc91f0078e9c1f82cfc1c841c..621ac275d9eb7fcf05763cf5c62d3436dda15e6b 100644
--- a/benchmarks/bench_translations.py
+++ b/benchmarks/bench_translations.py
@@ -1,13 +1,10 @@
 import numpy as np
 
-import pytest
-import pyopencl as cl
 from pyopencl.tools import (  # noqa
         pytest_generate_tests_for_pyopencl as pytest_generate_tests)
 
 from sumpy.expansion.multipole import (
         VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion,
-        VolumeTaylorMultipoleExpansionBase,
         LaplaceConformingVolumeTaylorMultipoleExpansion,
         HelmholtzConformingVolumeTaylorMultipoleExpansion)
 from sumpy.expansion.local import (
@@ -15,13 +12,11 @@ from sumpy.expansion.local import (
         LaplaceConformingVolumeTaylorLocalExpansion,
         HelmholtzConformingVolumeTaylorLocalExpansion)
 
-from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative,
-        DirectionalSourceDerivative)
+from sumpy.kernel import LaplaceKernel, HelmholtzKernel
 
 import logging
 logger = logging.getLogger(__name__)
 
-import sympy
 import six
 import pymbolic.mapper.flop_counter
 
@@ -29,6 +24,7 @@ import sumpy.symbolic as sym
 from sumpy.assignment_collection import SymbolicAssignmentCollection
 from sumpy.codegen import to_loopy_insns
 
+
 class Param:
     def __init__(self, dim, order):
         self.dim = dim
@@ -125,5 +121,3 @@ class Helmholtz2DTranslation(TranslationBenchmarkSuite):
         Param(2, 15),
         Param(2, 20),
     ]
-
-
diff --git a/contrib/translations/PDE-reduction-symbolic.ipynb b/contrib/translations/PDE-reduction-symbolic.ipynb
index bffcdf9e1e24beb09f0da9fbf50bf6d0ecdf1903..3550b0795bbf378f475db47e486472afd1e0ff48 100644
--- a/contrib/translations/PDE-reduction-symbolic.ipynb
+++ b/contrib/translations/PDE-reduction-symbolic.ipynb
@@ -125,7 +125,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.4"
+   "version": "3.6.5"
   }
  },
  "nbformat": 4,
diff --git a/sumpy/e2p.py b/sumpy/e2p.py
index 7b0072ad55708ba205ceed3448914ad0d8eda281..408b736408d9247db72d56da75e27e0b7dd6c857 100644
--- a/sumpy/e2p.py
+++ b/sumpy/e2p.py
@@ -91,7 +91,7 @@ class E2PBase(KernelCacheWrapper):
 
         coeff_exprs = [sym.Symbol("coeff%d" % i)
                 for i in range(len(self.expansion.get_coefficient_identifiers()))]
-        value = self.expansion.evaluate(coeff_exprs, bvec, rscale)
+        value = self.expansion.evaluate(coeff_exprs, bvec, rscale, sac=sac)
 
         result_names = [
             sac.assign_unique("result_%d_p" % i,
diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index 108bc67db749a50b17ac36a302675b4486df2787..85c93850795318acd8abda04755b93b343cddcbc 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -24,16 +24,16 @@ THE SOFTWARE.
 
 from six.moves import range
 import six
-import numpy as np
 import logging
 from pytools import memoize_method
 import sumpy.symbolic as sym
-from sumpy.tools import MiDerivativeTaker
 from collections import defaultdict
-
+from sumpy.tools import add_mi
+from .pde_utils import make_pde_sym, laplacian
 
 __doc__ = """
 .. autoclass:: ExpansionBase
+.. autoclass:: LinearPDEBasedExpansionTermsWrangler
 
 Expansion Factories
 ^^^^^^^^^^^^^^^^^^^
@@ -114,19 +114,21 @@ class ExpansionBase(object):
         """
         raise NotImplementedError
 
-    def coefficients_from_source(self, avec, bvec, rscale):
+    def coefficients_from_source(self, avec, bvec, rscale, sac=None):
         """Form an expansion from a source point.
 
         :arg avec: vector from source to center.
         :arg bvec: vector from center to target. Not usually necessary,
             except for line-Taylor expansion.
+        :arg sac: a symbolic assignment collection where temporary
+            expressions are stored.
 
         :returns: a list of :mod:`sympy` expressions representing
             the coefficients of the expansion.
         """
         raise NotImplementedError
 
-    def evaluate(self, coeffs, bvec, rscale):
+    def evaluate(self, coeffs, bvec, rscale, sac=None):
         """
         :return: a :mod:`sympy` expression corresponding
             to the evaluated expansion with the coefficients
@@ -158,26 +160,26 @@ class ExpansionBase(object):
 # }}}
 
 
-# {{{ derivative wrangler
+# {{{ expansion terms wrangler
+
+class ExpansionTermsWrangler(object):
 
-class DerivativeWrangler(object):
+    init_arg_names = ("order", "dim", "max_mi")
 
-    def __init__(self, order, dim):
+    def __init__(self, order, dim, max_mi=None):
         self.order = order
         self.dim = dim
+        self.max_mi = max_mi
 
     def get_coefficient_identifiers(self):
         raise NotImplementedError
 
     def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
-            rscale):
+            rscale, sac=None):
         raise NotImplementedError
 
     def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients,
-            rscale):
-        raise NotImplementedError
-
-    def get_derivative_taker(self, expr, var_list):
+            rscale, sac=None):
         raise NotImplementedError
 
     @memoize_method
@@ -190,19 +192,40 @@ class DerivativeWrangler(object):
                 as gnitstam)
 
         res = sorted(gnitstam(self.order, self.dim), key=sum)
+
+        def filter_tuple(tup):
+            if self.max_mi is None:
+                return True
+            for a, b in zip(tup, self.max_mi):
+                if a > b:
+                    return False
+            return True
+
+        res = list(filter(filter_tuple, res))
         return res
 
+    def copy(self, **kwargs):
+        new_kwargs = dict(
+                (name, getattr(self, name))
+                for name in self.init_arg_names)
+
+        for name in self.init_arg_names:
+            new_kwargs[name] = kwargs.pop(name, getattr(self, name))
 
-class FullDerivativeWrangler(DerivativeWrangler):
+        if kwargs:
+            raise TypeError("unexpected keyword arguments '%s'"
+                % ", ".join(kwargs))
 
-    def get_derivative_taker(self, expr, dvec):
-        return MiDerivativeTaker(expr, dvec)
+        return type(self)(**new_kwargs)
+
+
+class FullExpansionTermsWrangler(ExpansionTermsWrangler):
 
     get_coefficient_identifiers = (
-            DerivativeWrangler.get_full_coefficient_identifiers)
+            ExpansionTermsWrangler.get_full_coefficient_identifiers)
 
     def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
-            rscale):
+            rscale, sac=None):
         return stored_kernel_derivatives
 
     get_stored_mpole_coefficients_from_full = (
@@ -244,31 +267,80 @@ def _spmv(spmat, x, sparse_vectors):
 # }}}
 
 
-class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
+def _fast_spmv(reconstruct_matrix, vec, sac, transpose=False):
+    if not transpose:
+        res = [0] * len(reconstruct_matrix)
+        stored_idx = 0
+        for row, deps in enumerate(reconstruct_matrix):
+            if len(deps) == 0:
+                res[row] = vec[stored_idx]
+                stored_idx += 1
+            else:
+                for k, v in deps:
+                    res[row] += res[k] * v
+            new_sym = sym.Symbol(sac.assign_unique("expr", res[row]))
+            res[row] = new_sym
+        return res
+    else:
+        res = []
+        expr_all = list(vec)
+        for row, deps in reversed(list(enumerate(reconstruct_matrix))):
+            if len(deps) == 0:
+                res.append(expr_all[row])
+                continue
+            new_sym = sym.Symbol(sac.assign_unique("expr", expr_all[row]))
+            for k, v in deps:
+                expr_all[k] += new_sym * v
+        res.reverse()
+        return res
+
+
+class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler):
+    """
+    .. automethod:: __init__
+    .. automethod:: get_pde
+    .. automethod:: get_reduced_coeffs
+    """
+
+    init_arg_names = ("order", "dim", "max_mi")
+
+    def __init__(self, order, dim, max_mi=None):
+        r"""
+        :param order: order of the expansion
+        :param dim: number of dimensions
+        """
+        super(LinearPDEBasedExpansionTermsWrangler, self).__init__(order, dim,
+                max_mi)
 
     def get_coefficient_identifiers(self):
         return self.stored_identifiers
 
     def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
-            rscale):
-        coeff_matrix = self.get_coefficient_matrix(rscale)
-        return _spmv(coeff_matrix, stored_kernel_derivatives,
+            rscale, sac=None):
+        coeff_matrix, reconstruct_matrix, use_reconstruct = \
+            self.get_coefficient_matrix(rscale)
+        if not use_reconstruct or sac is None:
+            return _spmv(coeff_matrix, stored_kernel_derivatives,
                      sparse_vectors=False)
+        return _fast_spmv(reconstruct_matrix, stored_kernel_derivatives, sac)
 
     def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients,
-            rscale):
+            rscale, sac=None):
         # = M^T x, where M = coeff matrix
-
-        coeff_matrix = self.get_coefficient_matrix(rscale)
-        result = [0] * len(self.stored_identifiers)
-        for row, coeff in enumerate(full_mpole_coefficients):
-            for col, val in coeff_matrix[row]:
-                result[col] += coeff * val
-        return result
+        coeff_matrix, reconstruct_matrix, use_reconstruct = \
+            self.get_coefficient_matrix(rscale)
+        if not use_reconstruct or sac is None:
+            result = [0] * len(self.stored_identifiers)
+            for row, coeff in enumerate(full_mpole_coefficients):
+                for col, val in coeff_matrix[row]:
+                    result[col] += coeff * val
+            return result
+        return _fast_spmv(reconstruct_matrix, full_mpole_coefficients, sac,
+                transpose=True)
 
     @property
     def stored_identifiers(self):
-        stored_identifiers, coeff_matrix = self._get_stored_ids_and_coeff_mat()
+        stored_identifiers, coeff_matrix, _ = self.get_stored_ids_and_coeff_mat()
         return stored_identifiers
 
     @memoize_method
@@ -277,24 +349,26 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
         Return a matrix that expresses every derivative in terms of a
         set of "stored" derivatives.
 
-        For example, for the recurrence
+        For example, for the PDE::
 
             u_xx + u_yy + u_zz = 0
 
-        the coefficient matrix features the following entries:
+        the coefficient matrix features the following entries::
 
-              ... u_xx u_yy ... <= cols = only stored derivatives
-             ==================
-         ...| ...  ...  ... ...
-            |
-        u_zz| ...  -1   -1  ...
+                ... u_xx u_yy ... <= cols = only stored derivatives
+                ==================
+             ...| ...  ...  ... ...
+                |
+            u_zz| ...  -1   -1  ...
 
-           ^ rows = one for every derivative
+            ^ rows = one for every derivative
         """
-        stored_identifiers, coeff_matrix = self._get_stored_ids_and_coeff_mat()
+        stored_identifiers, coeff_matrix, reconstruct_matrix = \
+            self.get_stored_ids_and_coeff_mat()
 
         full_coeffs = self.get_full_coefficient_identifiers()
         matrix_rows = []
+        count_nonzero_coeff = -len(stored_identifiers)
         for irow, row in six.iteritems(coeff_matrix):
             # For eg: (u_xxx / rscale**3) = (u_yy / rscale**2) * coeff1 +
             #                               (u_xx / rscale**2) * coeff2
@@ -304,135 +378,159 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
             matrix_row = []
             for icol, coeff in row:
                 diff = row_rscale - sum(stored_identifiers[icol])
-                matrix_row.append((icol, coeff * rscale**diff))
+                mult = rscale**diff
+                matrix_row.append((icol, coeff * mult))
+            count_nonzero_coeff += len(row)
             matrix_rows.append((irow, matrix_row))
 
-        return defaultdict(list, matrix_rows)
+        if reconstruct_matrix is None:
+            return defaultdict(list, matrix_rows), None, False
 
-    @memoize_method
-    def _get_stored_ids_and_coeff_mat(self):
-        stored_identifiers = []
-        identifiers_so_far = {}
+        reconstruct_matrix_with_rscale = []
+        count_nonzero_reconstruct = 0
+        for row, deps in enumerate(reconstruct_matrix):
+            # For eg: (u_xxx / rscale**3) = (u_yy / rscale**2) * coeff1 +
+            #                               (u_xx / rscale**2) * coeff2
+            # is converted to u_xxx = u_yy * (rscale * coeff1) +
+            #                         u_xx * (rscale * coeff2)
+            row_rscale = sum(full_coeffs[row])
+            matrix_row = []
+            deps_with_rscale = []
+            for k, coeff in deps:
+                diff = row_rscale - sum(full_coeffs[k])
+                mult = rscale**diff
+                deps_with_rscale.append((k, coeff * mult))
+            count_nonzero_reconstruct += len(deps)
+            reconstruct_matrix_with_rscale.append(deps_with_rscale)
 
-        import time
-        start_time = time.time()
-        logger.debug("computing recurrence for Taylor coefficients: start")
+        use_reconstruct = count_nonzero_reconstruct < count_nonzero_coeff
 
-        # Sparse matrix, indexed by row
-        coeff_matrix_transpose = defaultdict(list)
+        return defaultdict(list, matrix_rows), reconstruct_matrix_with_rscale, \
+            use_reconstruct
 
-        # Build up the matrix transpose by row.
+    @memoize_method
+    def get_stored_ids_and_coeff_mat(self):
         from six import iteritems
-        for i, identifier in enumerate(self.get_full_coefficient_identifiers()):
-            expr = self.try_get_recurrence_for_derivative(
-                    identifier, identifiers_so_far)
-
-            if expr is None:
-                # Identifier should be stored
-                coeff_matrix_transpose[len(stored_identifiers)] = [(i, 1)]
-                stored_identifiers.append(identifier)
-            else:
-                expr = dict((identifiers_so_far[ident], val) for ident, val in
-                            iteritems(expr))
-                result = _spmv(coeff_matrix_transpose, expr, sparse_vectors=True)
-                for j, item in iteritems(result):
-                    coeff_matrix_transpose[j].append((i, item))
-
-            identifiers_so_far[identifier] = i
-
-        stored_identifiers = stored_identifiers
 
+        from pytools import ProcessLogger
+        plog = ProcessLogger(logger, "compute PDE for Taylor coefficients")
+
+        mis = self.get_full_coefficient_identifiers()
+        coeff_ident_enumerate_dict = dict((tuple(mi), i) for
+                                            (i, mi) in enumerate(mis))
+
+        pde_dict = self.get_pde().eq
+        for ident in pde_dict.keys():
+            if ident not in coeff_ident_enumerate_dict:
+                coeff_matrix = defaultdict(list)
+                reconstruct_matrix = []
+                for i in range(len(mis)):
+                    coeff_matrix[i] = [(i, 1)]
+                    reconstruct_matrix.append([])
+                return mis, coeff_matrix, reconstruct_matrix
+
+        max_mi_idx = max(coeff_ident_enumerate_dict[ident] for
+                         ident in pde_dict.keys())
+        max_mi = mis[max_mi_idx]
+        max_mi_coeff = pde_dict[max_mi]
+        max_mi_mult = -1/sym.sympify(max_mi_coeff)
+
+        def is_stored(mi):
+            """
+            A multi_index mi is not stored if mi >= max_mi
+            """
+            return any(mi[d] < max_mi[d] for d in range(self.dim))
+
+        stored_identifiers = [mi for mi in mis if is_stored(mi)]
+        stored_ident_enumerate_dict = dict((tuple(mi), i) for
+                                            (i, mi) in enumerate(stored_identifiers))
+
+        coeff_matrix_dict = defaultdict(lambda: defaultdict(lambda: 0))
+        reconstruct_matrix = []
+        for i, mi in enumerate(mis):
+            reconstruct_matrix.append([])
+            if is_stored(mi):
+                coeff_matrix_dict[i][stored_ident_enumerate_dict[mi]] = 1
+                continue
+            diff = [mi[d] - max_mi[d] for d in range(self.dim)]
+            for other_mi, coeff in iteritems(pde_dict):
+                j = coeff_ident_enumerate_dict[add_mi(other_mi, diff)]
+                if i == j:
+                    continue
+                # PDE might not have max_mi_coeff = -1, divide by -max_mi_coeff
+                # to get a relation of the form, u_zz = - u_xx - u_yy for Laplace 3D.
+                reconstruct_matrix[i].append((j, coeff*max_mi_mult))
+                # j can be a stored or a non-stored multi-index
+                # Look at the coeff_matrix of j to get the j as a linear combination
+                # of stored coefficients.
+                for dep, other_coeff in iteritems(coeff_matrix_dict[j]):
+                    coeff_matrix_dict[i][dep] += other_coeff*coeff*max_mi_mult
+
+        # coeff_matrix is a dictionary of lists. Each key in the dictionary
+        # is a row number and the values are pairs of column number and value.
         coeff_matrix = defaultdict(list)
-        for i, row in iteritems(coeff_matrix_transpose):
-            for j, val in row:
-                coeff_matrix[j].append((i, val))
+        for row, deps in iteritems(coeff_matrix_dict):
+            for col, val in iteritems(deps):
+                if val != 0:
+                    coeff_matrix[row].append((col, val))
 
-        logger.debug("computing recurrence for Taylor coefficients: "
-                     "done after {dur:.2f} seconds"
-                     .format(dur=time.time() - start_time))
+        plog.done()
 
-        logger.debug("number of Taylor coefficients was reduced from {orig} to {red}"
+        print("number of Taylor coefficients was reduced from {orig} to {red}"
                      .format(orig=len(self.get_full_coefficient_identifiers()),
                              red=len(stored_identifiers)))
 
-        return stored_identifiers, coeff_matrix
+        return stored_identifiers, coeff_matrix, reconstruct_matrix
 
-    def try_get_recurrence_for_derivative(self, deriv, in_terms_of):
-        """
-        :arg deriv: a tuple of integers identifying a derivative for which
-            a recurrence is sought
-        :arg in_terms_of: a container supporting efficient containment testing
-            indicating availability of derivatives for use in recurrences
+    def get_pde(self):
+        r"""
+        Returns a PDE. A PDE stores a dictionary of (mi, coeff)
+        where mi is the multi-index of the  derivative and coeff is the
+        coefficient
         """
 
         raise NotImplementedError
 
-    def get_derivative_taker(self, expr, var_list):
-        from sumpy.tools import LinearRecurrenceBasedMiDerivativeTaker
-        return LinearRecurrenceBasedMiDerivativeTaker(expr, var_list, self)
 
+class LaplaceExpansionTermsWrangler(LinearPDEBasedExpansionTermsWrangler):
 
-class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
+    init_arg_names = ("order", "dim", "max_mi")
 
-    def try_get_recurrence_for_derivative(self, deriv, in_terms_of):
-        deriv = np.array(deriv, dtype=int)
+    def __init__(self, order, dim, max_mi=None):
+        super(LaplaceExpansionTermsWrangler, self).__init__(order=order, dim=dim,
+            max_mi=max_mi)
 
-        for dim in np.where(2 <= deriv)[0]:
-            # Check if we can reduce this dimension in terms of the other
-            # dimensions.
+    def get_pde(self):
+        w = make_pde_sym(self.dim)
+        return laplacian(w)
 
-            reduced_deriv = deriv.copy()
-            reduced_deriv[dim] -= 2
-            coeffs = {}
-
-            for other_dim in range(self.dim):
-                if other_dim == dim:
-                    continue
-                needed_deriv = reduced_deriv.copy()
-                needed_deriv[other_dim] += 2
-                needed_deriv = tuple(needed_deriv)
-
-                if needed_deriv not in in_terms_of:
-                    break
-
-                coeffs[needed_deriv] = -1
-            else:
-                return coeffs
 
+class HelmholtzExpansionTermsWrangler(LinearPDEBasedExpansionTermsWrangler):
 
-class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
+    init_arg_names = ("order", "dim", "helmholtz_k_name", "max_mi")
 
-    def __init__(self, order, dim, helmholtz_k_name):
-        super(HelmholtzDerivativeWrangler, self).__init__(order, dim)
+    def __init__(self, order, dim, helmholtz_k_name, max_mi=None):
         self.helmholtz_k_name = helmholtz_k_name
+        super(HelmholtzExpansionTermsWrangler, self).__init__(order=order, dim=dim,
+            max_mi=max_mi)
 
-    def try_get_recurrence_for_derivative(self, deriv, in_terms_of):
-        deriv = np.array(deriv, dtype=int)
+    def get_pde(self, **kwargs):
+        w = make_pde_sym(self.dim)
+        k = sym.Symbol(self.helmholtz_k_name)
+        return (laplacian(w) + k**2 * w)
 
-        for dim in np.where(2 <= deriv)[0]:
-            # Check if we can reduce this dimension in terms of the other
-            # dimensions.
 
-            reduced_deriv = deriv.copy()
-            reduced_deriv[dim] -= 2
-            coeffs = {}
+class BiharmonicExpansionTermsWrangler(LinearPDEBasedExpansionTermsWrangler):
 
-            for other_dim in range(self.dim):
-                if other_dim == dim:
-                    continue
-                needed_deriv = reduced_deriv.copy()
-                needed_deriv[other_dim] += 2
-                needed_deriv = tuple(needed_deriv)
+    init_arg_names = ("order", "dim", "max_mi")
 
-                if needed_deriv not in in_terms_of:
-                    break
-
-                coeffs[needed_deriv] = -1
-            else:
-                k = sym.Symbol(self.helmholtz_k_name)
-                coeffs[tuple(reduced_deriv)] = -k*k
-                return coeffs
+    def __init__(self, order, dim, max_mi=None):
+        super(BiharmonicExpansionTermsWrangler, self).__init__(order=order, dim=dim,
+            max_mi=max_mi)
 
+    def get_pde(self, **kwargs):
+        w = make_pde_sym(self.dim)
+        return laplacian(laplacian(w))
 # }}}
 
 
@@ -441,31 +539,33 @@ class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
 class VolumeTaylorExpansionBase(object):
 
     @classmethod
-    def get_or_make_derivative_wrangler(cls, *key):
+    def get_or_make_expansion_terms_wrangler(cls, *key):
         """
-        This stores the derivative wrangler at the class attribute level because
-        precomputing the derivative wrangler may become expensive.
+        This stores the expansion terms wrangler at the class attribute level
+        because recreating the expansion terms wrangler implicitly empties its
+        caches.
         """
         try:
-            wrangler = cls.derivative_wrangler_cache[key]
+            wrangler = cls.expansion_terms_wrangler_cache[key]
         except KeyError:
-            wrangler = cls.derivative_wrangler_class(*key)
-            cls.derivative_wrangler_cache[key] = wrangler
+            wrangler = cls.expansion_terms_wrangler_class(*key)
+            cls.expansion_terms_wrangler_cache[key] = wrangler
 
         return wrangler
 
     @property
-    def derivative_wrangler(self):
-        return self.get_or_make_derivative_wrangler(*self.derivative_wrangler_key)
+    def expansion_terms_wrangler(self):
+        return self.get_or_make_expansion_terms_wrangler(
+                *self.expansion_terms_wrangler_key)
 
     def get_coefficient_identifiers(self):
         """
         Returns the identifiers of the coefficients that actually get stored.
         """
-        return self.derivative_wrangler.get_coefficient_identifiers()
+        return self.expansion_terms_wrangler.get_coefficient_identifiers()
 
     def get_full_coefficient_identifiers(self):
-        return self.derivative_wrangler.get_full_coefficient_identifiers()
+        return self.expansion_terms_wrangler.get_full_coefficient_identifiers()
 
     @property
     @memoize_method
@@ -479,33 +579,43 @@ class VolumeTaylorExpansionBase(object):
 
 class VolumeTaylorExpansion(VolumeTaylorExpansionBase):
 
-    derivative_wrangler_class = FullDerivativeWrangler
-    derivative_wrangler_cache = {}
+    expansion_terms_wrangler_class = FullExpansionTermsWrangler
+    expansion_terms_wrangler_cache = {}
 
     # not user-facing, be strict about having to pass use_rscale
     def __init__(self, kernel, order, use_rscale):
-        self.derivative_wrangler_key = (order, kernel.dim)
+        self.expansion_terms_wrangler_key = (order, kernel.dim)
 
 
 class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
 
-    derivative_wrangler_class = LaplaceDerivativeWrangler
-    derivative_wrangler_cache = {}
+    expansion_terms_wrangler_class = LaplaceExpansionTermsWrangler
+    expansion_terms_wrangler_cache = {}
 
     # not user-facing, be strict about having to pass use_rscale
     def __init__(self, kernel, order, use_rscale):
-        self.derivative_wrangler_key = (order, kernel.dim)
+        self.expansion_terms_wrangler_key = (order, kernel.dim)
 
 
 class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
 
-    derivative_wrangler_class = HelmholtzDerivativeWrangler
-    derivative_wrangler_cache = {}
+    expansion_terms_wrangler_class = HelmholtzExpansionTermsWrangler
+    expansion_terms_wrangler_cache = {}
 
     # not user-facing, be strict about having to pass use_rscale
     def __init__(self, kernel, order, use_rscale):
         helmholtz_k_name = kernel.get_base_kernel().helmholtz_k_name
-        self.derivative_wrangler_key = (order, kernel.dim, helmholtz_k_name)
+        self.expansion_terms_wrangler_key = (order, kernel.dim, helmholtz_k_name)
+
+
+class BiharmonicConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
+
+    expansion_terms_wrangler_class = BiharmonicExpansionTermsWrangler
+    expansion_terms_wrangler_cache = {}
+
+    # not user-facing, be strict about having to pass use_rscale
+    def __init__(self, kernel, order, use_rscale):
+        self.expansion_terms_wrangler_key = (order, kernel.dim)
 
 # }}}
 
@@ -555,49 +665,58 @@ class DefaultExpansionFactory(ExpansionFactoryBase):
     def get_local_expansion_class(self, base_kernel):
         """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*.
         """
-        from sumpy.kernel import HelmholtzKernel, LaplaceKernel, YukawaKernel
+        from sumpy.kernel import (HelmholtzKernel, LaplaceKernel, YukawaKernel,
+                BiharmonicKernel, StokesletKernel, StressletKernel)
+
+        from sumpy.expansion.local import (H2DLocalExpansion, Y2DLocalExpansion,
+                HelmholtzConformingVolumeTaylorLocalExpansion,
+                LaplaceConformingVolumeTaylorLocalExpansion,
+                BiharmonicConformingVolumeTaylorLocalExpansion,
+                VolumeTaylorLocalExpansion)
+
         if (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel)
                 and base_kernel.dim == 2):
-            from sumpy.expansion.local import H2DLocalExpansion
             return H2DLocalExpansion
         elif (isinstance(base_kernel.get_base_kernel(), YukawaKernel)
                 and base_kernel.dim == 2):
-            from sumpy.expansion.local import Y2DLocalExpansion
             return Y2DLocalExpansion
         elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel):
-            from sumpy.expansion.local import \
-                    HelmholtzConformingVolumeTaylorLocalExpansion
             return HelmholtzConformingVolumeTaylorLocalExpansion
         elif isinstance(base_kernel.get_base_kernel(), LaplaceKernel):
-            from sumpy.expansion.local import \
-                    LaplaceConformingVolumeTaylorLocalExpansion
             return LaplaceConformingVolumeTaylorLocalExpansion
+        elif isinstance(base_kernel.get_base_kernel(),
+                (BiharmonicKernel, StokesletKernel, StressletKernel)):
+            return BiharmonicConformingVolumeTaylorLocalExpansion
         else:
-            from sumpy.expansion.local import VolumeTaylorLocalExpansion
             return VolumeTaylorLocalExpansion
 
     def get_multipole_expansion_class(self, base_kernel):
         """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*.
         """
-        from sumpy.kernel import HelmholtzKernel, LaplaceKernel, YukawaKernel
+        from sumpy.kernel import (HelmholtzKernel, LaplaceKernel, YukawaKernel,
+                BiharmonicKernel, StokesletKernel, StressletKernel)
+
+        from sumpy.expansion.multipole import (H2DMultipoleExpansion,
+                Y2DMultipoleExpansion,
+                LaplaceConformingVolumeTaylorMultipoleExpansion,
+                HelmholtzConformingVolumeTaylorMultipoleExpansion,
+                BiharmonicConformingVolumeTaylorMultipoleExpansion,
+                VolumeTaylorMultipoleExpansion)
+
         if (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel)
                 and base_kernel.dim == 2):
-            from sumpy.expansion.multipole import H2DMultipoleExpansion
             return H2DMultipoleExpansion
         elif (isinstance(base_kernel.get_base_kernel(), YukawaKernel)
                 and base_kernel.dim == 2):
-            from sumpy.expansion.multipole import Y2DMultipoleExpansion
             return Y2DMultipoleExpansion
         elif isinstance(base_kernel.get_base_kernel(), LaplaceKernel):
-            from sumpy.expansion.multipole import (
-                    LaplaceConformingVolumeTaylorMultipoleExpansion)
             return LaplaceConformingVolumeTaylorMultipoleExpansion
         elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel):
-            from sumpy.expansion.multipole import (
-                    HelmholtzConformingVolumeTaylorMultipoleExpansion)
             return HelmholtzConformingVolumeTaylorMultipoleExpansion
+        elif isinstance(base_kernel.get_base_kernel(),
+                (BiharmonicKernel, StokesletKernel, StressletKernel)):
+            return BiharmonicConformingVolumeTaylorMultipoleExpansion
         else:
-            from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion
             return VolumeTaylorMultipoleExpansion
 
 # }}}
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index b4b435de5ce1aa778b5c176b9a6df18b3d88c8dc..a58c598173a352afaa92d6a6d06855ca07275457 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -27,7 +27,8 @@ import sumpy.symbolic as sym
 
 from sumpy.expansion import (
     ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion,
-    HelmholtzConformingVolumeTaylorExpansion)
+    HelmholtzConformingVolumeTaylorExpansion,
+    BiharmonicConformingVolumeTaylorExpansion)
 
 
 class LocalExpansionBase(ExpansionBase):
@@ -57,7 +58,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase):
     def get_coefficient_identifiers(self):
         return list(range(self.order+1))
 
-    def coefficients_from_source(self, avec, bvec, rscale):
+    def coefficients_from_source(self, avec, bvec, rscale, sac=None):
         # no point in heeding rscale here--just ignore it
         if bvec is None:
             raise RuntimeError("cannot use line-Taylor expansions in a setting "
@@ -98,7 +99,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase):
                     .subs("tau", 0)
                     for i in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec, rscale):
+    def evaluate(self, coeffs, bvec, rscale, sac=None):
         # no point in heeding rscale here--just ignore it
         from pytools import factorial
         return sym.Add(*(
@@ -115,7 +116,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
     Coefficients represent derivative values of the kernel.
     """
 
-    def coefficients_from_source(self, avec, bvec, rscale):
+    def coefficients_from_source(self, avec, bvec, rscale, sac=None):
         from sumpy.tools import MiDerivativeTaker
         ppkernel = self.kernel.postprocess_at_source(
                 self.kernel.get_expression(avec), avec)
@@ -125,11 +126,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
                 taker.diff(mi) * rscale ** sum(mi)
                 for mi in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec, rscale):
+    def evaluate(self, coeffs, bvec, rscale, sac=None):
         from pytools import factorial
 
         evaluated_coeffs = (
-            self.derivative_wrangler.get_full_kernel_derivatives_from_stored(
+            self.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored(
                 coeffs, rscale))
 
         bvec = [b*rscale**-1 for b in bvec]
@@ -239,23 +240,62 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
 
             from sumpy.tools import add_mi
 
+            max_mi = [0]*self.dim
+            for i in range(self.dim):
+                max_mi[i] = max(mi[i] for mi in
+                                  src_expansion.get_coefficient_identifiers())
+                max_mi[i] += max(mi[i] for mi in
+                                  self.get_coefficient_identifiers())
+
+            # Create a expansion terms wrangler for derivatives up to order
+            # (tgt order)+(src order) including a corresponding reduction matrix
+            tgtplusderiv_exp_terms_wrangler = \
+                src_expansion.expansion_terms_wrangler.copy(
+                        order=self.order + src_expansion.order, max_mi=tuple(max_mi))
+            tgtplusderiv_coeff_ids = \
+                tgtplusderiv_exp_terms_wrangler.get_coefficient_identifiers()
+            tgtplusderiv_full_coeff_ids = \
+                tgtplusderiv_exp_terms_wrangler.get_full_coefficient_identifiers()
+
+            tgtplusderiv_ident_to_index = dict((ident, i) for i, ident in
+                                enumerate(tgtplusderiv_full_coeff_ids))
+
             result = []
-            for deriv in self.get_coefficient_identifiers():
-                local_result = []
+            for lexp_mi in self.get_coefficient_identifiers():
+                lexp_mi_terms = []
+
+                # Embed the source coefficients in the coefficient array
+                # for the (tgt order)+(src order) wrangler, offset by lexp_mi.
+                embedded_coeffs = [0] * len(tgtplusderiv_full_coeff_ids)
                 for coeff, term in zip(
                         src_coeff_exprs,
                         src_expansion.get_coefficient_identifiers()):
-
+                    embedded_coeffs[
+                            tgtplusderiv_ident_to_index[add_mi(lexp_mi, term)]] \
+                                    = coeff
+
+                # Compress the embedded coefficient set
+                stored_coeffs = tgtplusderiv_exp_terms_wrangler \
+                        .get_stored_mpole_coefficients_from_full(
+                                embedded_coeffs, src_rscale)
+
+                # Sum the embedded coefficient set
+                for i, coeff in enumerate(stored_coeffs):
+                    if coeff == 0:
+                        continue
+                    nderivatives_for_scaling = \
+                            sum(tgtplusderiv_coeff_ids[i])-sum(lexp_mi)
                     kernel_deriv = (
                             src_expansion.get_scaled_multipole(
-                                taker.diff(add_mi(deriv, term)),
+                                taker.diff(tgtplusderiv_coeff_ids[i]),
                                 dvec, src_rscale,
-                                nderivatives=sum(deriv) + sum(term),
-                                nderivatives_for_scaling=sum(term)))
+                                nderivatives=sum(tgtplusderiv_coeff_ids[i]),
+                                nderivatives_for_scaling=nderivatives_for_scaling))
+
+                    lexp_mi_terms.append(
+                            coeff * kernel_deriv * tgt_rscale**sum(lexp_mi))
+                result.append(sym.Add(*lexp_mi_terms))
 
-                    local_result.append(
-                            coeff * kernel_deriv * tgt_rscale**sum(deriv))
-                result.append(sym.Add(*local_result))
         else:
             from sumpy.tools import MiDerivativeTaker
             # Rscale/operand magnitude is fairly sensitive to the order of
@@ -307,6 +347,16 @@ class HelmholtzConformingVolumeTaylorLocalExpansion(
         HelmholtzConformingVolumeTaylorExpansion.__init__(
                 self, kernel, order, use_rscale)
 
+
+class BiharmonicConformingVolumeTaylorLocalExpansion(
+        BiharmonicConformingVolumeTaylorExpansion,
+        VolumeTaylorLocalExpansionBase):
+
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale)
+        BiharmonicConformingVolumeTaylorExpansion.__init__(
+                self, kernel, order, use_rscale)
+
 # }}}
 
 
@@ -319,7 +369,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
     def get_coefficient_identifiers(self):
         return list(range(-self.order, self.order+1))
 
-    def coefficients_from_source(self, avec, bvec, rscale):
+    def coefficients_from_source(self, avec, bvec, rscale, sac=None):
         if not self.use_rscale:
             rscale = 1
 
@@ -337,7 +387,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
                     * sym.exp(sym.I * l * source_angle_rel_center), avec)
                     for l in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec, rscale):
+    def evaluate(self, coeffs, bvec, rscale, sac=None):
         if not self.use_rscale:
             rscale = 1
 
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 068fea2393875f7b58904f088f7d58e8939e09b3..1a53aa40f15f78fb3d94c1b2f8bb9932d6daba90 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -28,7 +28,8 @@ import sumpy.symbolic as sym  # noqa
 from sumpy.symbolic import vector_xreplace
 from sumpy.expansion import (
     ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion,
-    HelmholtzConformingVolumeTaylorExpansion)
+    HelmholtzConformingVolumeTaylorExpansion,
+    BiharmonicConformingVolumeTaylorExpansion)
 from pytools import cartesian_product
 
 import logging
@@ -55,7 +56,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
     Coefficients represent the terms in front of the kernel derivatives.
     """
 
-    def coefficients_from_source(self, avec, bvec, rscale):
+    def coefficients_from_source(self, avec, bvec, rscale, sac=None):
         from sumpy.kernel import DirectionalSourceDerivative
         kernel = self.kernel
 
@@ -107,8 +108,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
                     mi_power(avec, mi) / mi_factorial(mi)
                     for mi in self.get_full_coefficient_identifiers()]
         return (
-            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
-                result, rscale))
+            self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full(
+                result, rscale, sac=sac))
 
     def get_scaled_multipole(self, expr, bvec, rscale, nderivatives,
             nderivatives_for_scaling=None):
@@ -126,7 +127,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
         else:
             return (rscale**nderivatives_for_scaling * expr)
 
-    def evaluate(self, coeffs, bvec, rscale):
+    def evaluate(self, coeffs, bvec, rscale, sac=None):
         if not self.use_rscale:
             rscale = 1
 
@@ -140,8 +141,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
         return result
 
     def get_kernel_derivative_taker(self, bvec):
-        return (self.derivative_wrangler.get_derivative_taker(
-            self.kernel.get_expression(bvec), bvec))
+        from sumpy.tools import MiDerivativeTaker
+        return MiDerivativeTaker(self.kernel.get_expression(bvec), bvec)
 
     def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
             dvec, tgt_rscale):
@@ -239,7 +240,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
 
         logger.info("building translation operator: done")
         return (
-            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
+            self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full(
                 result, tgt_rscale))
 
 
@@ -271,6 +272,16 @@ class HelmholtzConformingVolumeTaylorMultipoleExpansion(
         HelmholtzConformingVolumeTaylorExpansion.__init__(
                 self, kernel, order, use_rscale)
 
+
+class BiharmonicConformingVolumeTaylorMultipoleExpansion(
+        BiharmonicConformingVolumeTaylorExpansion,
+        VolumeTaylorMultipoleExpansionBase):
+
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale)
+        BiharmonicConformingVolumeTaylorExpansion.__init__(
+                self, kernel, order, use_rscale)
+
 # }}}
 
 
@@ -283,7 +294,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase):
     def get_coefficient_identifiers(self):
         return list(range(-self.order, self.order+1))
 
-    def coefficients_from_source(self, avec, bvec, rscale):
+    def coefficients_from_source(self, avec, bvec, rscale, sac=None):
         if not self.use_rscale:
             rscale = 1
 
@@ -303,7 +314,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase):
                     avec)
                 for l in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec, rscale):
+    def evaluate(self, coeffs, bvec, rscale, sac=None):
         if not self.use_rscale:
             rscale = 1
 
diff --git a/sumpy/expansion/pde_utils.py b/sumpy/expansion/pde_utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..146bd510183909adf8ca704ac269fc4a95a4718d
--- /dev/null
+++ b/sumpy/expansion/pde_utils.py
@@ -0,0 +1,91 @@
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2019 Isuru Fernando"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+from sumpy.tools import add_mi
+
+
+class PDE(object):
+    r"""
+    Represents a iscalar PDEs of dimension `dim`. It is represented by a
+    dictionary. A dictionary maps a a multi-index given as a tuple
+    to the coefficient.
+    """
+    def __init__(self, dim, eq):
+        """
+        :arg dim: dimension of the PDE
+        :arg eq: A dictionary mapping a multi-index to a value
+        """
+        self.dim = dim
+        self.eq = eq
+
+    def __mul__(self, param):
+        res = PDE(self.dim, {})
+        for k, v in self.eq.items():
+            res.eq[k] = v * param
+        return res
+
+    __rmul__ = __mul__
+
+    def __add__(self, other_pde):
+        assert self.dim == other_pde.dim
+        res = PDE(self.dim, self.eq.copy())
+        for k, v in other_pde.eq.items():
+            if k in res.eq:
+                res.eq[k] += v
+            else:
+                res.eq[k] = v
+        return res
+
+    __radd__ = __add__
+
+    def __sub__(self, other_pde):
+        return self + (-1)*other_pde
+
+    def __repr__(self):
+        return repr(self.eq)
+
+
+def laplacian(pde):
+    dim = pde.dim
+    res = PDE(dim, {})
+    for j in range(dim):
+        mi = [0]*dim
+        mi[j] = 2
+        res = res + diff(pde, tuple(mi))
+    return res
+
+
+def diff(pde, mi):
+    res = PDE(pde.dim, {})
+    for eq_mi, v in pde.eq.items():
+        res.eq[add_mi(eq_mi, mi)] = v
+    return res
+
+
+def make_pde_sym(dim):
+    """
+    Returns a PDE u = 0
+    """
+    mi = tuple([0]*dim)
+    return PDE(dim, {mi: 1})
diff --git a/sumpy/p2e.py b/sumpy/p2e.py
index daa7d93dc5b393d358424031500a4915e2c2dd9f..f5472fc42e306c21eeec5091b26d548cdb66c71a 100644
--- a/sumpy/p2e.py
+++ b/sumpy/p2e.py
@@ -86,9 +86,9 @@ class P2EBase(KernelCacheWrapper):
         sac = SymbolicAssignmentCollection()
 
         coeff_names = [
-                sac.assign_unique("coeff%d" % i, coeff_i)
-                for i, coeff_i in enumerate(
-                    self.expansion.coefficients_from_source(avec, None, rscale))]
+            sac.assign_unique("coeff%d" % i, coeff_i)
+            for i, coeff_i in enumerate(
+                self.expansion.coefficients_from_source(avec, None, rscale, sac))]
 
         sac.run_global_cse()
 
diff --git a/sumpy/tools.py b/sumpy/tools.py
index 4d1098429d1c9db43af4a6e26fbb63509438a375..2719f43a2e128da47997f0d1a010072d9bcec1dd 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -106,53 +106,6 @@ class MiDerivativeTaker(object):
                 if (np.array(mi) >= np.array(other_mi)).all()),
             key=lambda other_mi: sum(self.mi_dist(mi, other_mi)))
 
-
-class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker):
-    """
-    The derivative taker for expansions that use
-    :class:`sumpy.expansion.LinearRecurrenceBasedDerivativeWrangler`
-    """
-
-    def __init__(self, expr, var_list, wrangler):
-        super(LinearRecurrenceBasedMiDerivativeTaker, self).__init__(
-                expr, var_list)
-        self.wrangler = wrangler
-
-    @memoize_method
-    def diff(self, mi):
-        """
-        :arg mi: a multi-index (tuple) indicating how many x/y derivatives are
-            to be taken.
-        """
-        try:
-            expr = self.cache_by_mi[mi]
-        except KeyError:
-            from six import iteritems
-            from sumpy.symbolic import Add
-
-            closest_mi = self.get_closest_cached_mi(mi)
-            expr = self.cache_by_mi[closest_mi]
-
-            # Try to reduce the derivative using recurrences first, and if that
-            # fails fall back to derivative taking.
-            for next_deriv, next_mi in (
-                        self.get_derivative_taking_sequence(closest_mi, mi)):
-
-                recurrence = (
-                        self.wrangler.try_get_recurrence_for_derivative(
-                            next_mi, self.cache_by_mi))
-
-                if recurrence is not None:
-                    expr = Add(*tuple(
-                            coeff * self.cache_by_mi[ident]
-                            for ident, coeff in iteritems(recurrence)))
-                else:
-                    expr = expr.diff(next_deriv)
-
-                self.cache_by_mi[next_mi] = expr
-
-        return expr
-
 # }}}
 
 
@@ -717,4 +670,94 @@ def my_syntactic_subs(expr, subst_dict):
         return expr
 
 
+def reduced_row_echelon_form(m):
+    """Calculates a reduced row echelon form of a
+    matrix `m`.
+
+    :arg m: a 2D :class:`numpy.ndarray` or a list of lists or a sympy Matrix
+    :return: reduced row echelon form as a 2D :class:`numpy.ndarray`
+             and a list of pivots
+    """
+
+    mat = np.array(m, dtype=object)
+    index = 0
+    nrows = mat.shape[0]
+    ncols = mat.shape[1]
+    pivot_cols = []
+    for i in range(ncols):
+        if index == nrows:
+            break
+        pivot = nrows
+        for k in range(index, nrows):
+            if mat[k, i] != 0 and pivot == nrows:
+                pivot = k
+            if abs(mat[k, i]) == 1:
+                pivot = k
+                break
+        if pivot == nrows:
+            continue
+        if pivot != index:
+            mat[[pivot, index], :] = mat[[index, pivot], :]
+
+        pivot_cols.append(i)
+        scale = mat[index, i]
+        if isinstance(scale, (int, sym.Integer)):
+            scale = int(scale)
+
+        for j in range(mat.shape[1]):
+            elem = mat[index, j]
+            if isinstance(scale, int) and isinstance(elem, (int, sym.Integer)):
+                quo = int(elem) // scale
+                if quo * scale == elem:
+                    mat[index, j] = quo
+                    continue
+            mat[index, j] = sym.sympify(elem)/scale
+
+        for j in range(nrows):
+            if (j == index):
+                continue
+
+            scale = mat[j, i]
+            if scale != 0:
+                mat[j, :] = mat[j, :] - mat[index, :]*scale
+
+        index = index + 1
+
+    return mat, pivot_cols
+
+
+def nullspace(m):
+    """Calculates the nullspace of a matrix `m`.
+
+    :arg m: a 2D :class:`numpy.ndarray` or a list of lists or a sympy Matrix
+    :return: nullspace of `m` as a 2D :class:`numpy.ndarray`
+    """
+    mat, pivot_cols = reduced_row_echelon_form(m)
+    pivot_cols = list(pivot_cols)
+    cols = mat.shape[1]
+
+    free_vars = [i for i in range(cols) if i not in pivot_cols]
+
+    n = []
+    for free_var in free_vars:
+        vec = [0]*cols
+        vec[free_var] = 1
+        for piv_row, piv_col in enumerate(pivot_cols):
+            for pos in pivot_cols[piv_row+1:] + [free_var]:
+                if isinstance(mat[piv_row, pos], sym.Integer):
+                    vec[piv_col] -= int(mat[piv_row, pos])
+                else:
+                    vec[piv_col] -= mat[piv_row, pos]
+        n.append(vec)
+    return np.array(n, dtype=object).T
+
+
+def solve_symbolic(A, b):  # noqa: N803
+    if isinstance(A, sym.Matrix):
+        big = A.row_join(b)
+    else:
+        big = np.hstack((A, b))
+    red = reduced_row_echelon_form(big)[0]
+    return red[:, big.shape[0]:]
+
 # vim: fdm=marker
diff --git a/test/test_kernels.py b/test/test_kernels.py
index 834d823321bec0c05b3253e4c8718ff9145b343f..8225723cdd4f7df2afb7c487a578cc8c9ee81d28 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -22,8 +22,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-from six.moves import range
-
 import numpy as np
 import numpy.linalg as la
 import sys
@@ -37,13 +35,15 @@ from sumpy.expansion.multipole import (
         VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion,
         VolumeTaylorMultipoleExpansionBase,
         LaplaceConformingVolumeTaylorMultipoleExpansion,
-        HelmholtzConformingVolumeTaylorMultipoleExpansion)
+        HelmholtzConformingVolumeTaylorMultipoleExpansion,
+        BiharmonicConformingVolumeTaylorMultipoleExpansion)
 from sumpy.expansion.local import (
         VolumeTaylorLocalExpansion, H2DLocalExpansion,
         LaplaceConformingVolumeTaylorLocalExpansion,
-        HelmholtzConformingVolumeTaylorLocalExpansion)
+        HelmholtzConformingVolumeTaylorLocalExpansion,
+        BiharmonicConformingVolumeTaylorLocalExpansion)
 from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative,
-        DirectionalSourceDerivative, BiharmonicKernel)
+        DirectionalSourceDerivative, BiharmonicKernel, StokesletKernel)
 from pytools.convergence import PConvergenceVerifier
 
 import logging
@@ -160,6 +160,8 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative):
             extra_kwargs["k"] = 0.2 * (0.707 + 0.707j)
         else:
             extra_kwargs["k"] = 0.2
+    if isinstance(base_knl, StokesletKernel):
+        extra_kwargs["mu"] = 0.2
 
     if with_source_derivative:
         knl = DirectionalSourceDerivative(base_knl, "dir_vec")
@@ -187,7 +189,7 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative):
     else:
         h_values = [1/2, 1/3, 1/5]
 
-    center = np.array([2, 1], np.float64)
+    center = np.array([2, 1, 0][:knl.dim], np.float64)
     sources = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64))
             + center[:, np.newaxis])
 
@@ -207,15 +209,15 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative):
 
     for h in h_values:
         if issubclass(expn_class, LocalExpansionBase):
-            loc_center = np.array([5.5, 0.0]) + center
+            loc_center = np.array([5.5, 0.0, 0.0][:knl.dim]) + center
             centers = np.array(loc_center, dtype=np.float64).reshape(knl.dim, 1)
             fp = FieldPlotter(loc_center, extent=h, npoints=res)
         else:
-            eval_center = np.array([1/h, 0.0]) + center
+            eval_center = np.array([1/h, 0.0, 0.0][:knl.dim]) + center
             fp = FieldPlotter(eval_center, extent=0.1, npoints=res)
-            centers = (
-                    np.array([0.0, 0.0], dtype=np.float64).reshape(knl.dim, 1)
-                    + center[:, np.newaxis])
+            centers = (np.array([0.0, 0.0, 0.0][:knl.dim],
+                                dtype=np.float64).reshape(knl.dim, 1)
+                        + center[:, np.newaxis])
 
         targets = fp.points
 
@@ -351,7 +353,11 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative):
     (HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
     (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorLocalExpansion,
      HelmholtzConformingVolumeTaylorMultipoleExpansion),
-    (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion)
+    (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion),
+    (StokesletKernel(2, 0, 0), VolumeTaylorLocalExpansion,
+     VolumeTaylorMultipoleExpansion),
+    (StokesletKernel(2, 0, 0), BiharmonicConformingVolumeTaylorLocalExpansion,
+     BiharmonicConformingVolumeTaylorMultipoleExpansion),
     ])
 def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
     logging.basicConfig(level=logging.INFO)
@@ -372,9 +378,11 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
     extra_kwargs = {}
     if isinstance(knl, HelmholtzKernel):
         extra_kwargs["k"] = 0.05
+    if isinstance(knl, StokesletKernel):
+        extra_kwargs["mu"] = 0.05
 
     # Just to make sure things also work away from the origin
-    origin = np.array([2, 1], np.float64)
+    origin = np.array([2, 1, 0][:knl.dim], np.float64)
     sources = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64))
             + origin[:, np.newaxis])
     strengths = np.ones(nsources, dtype=np.float64) * (1/nsources)
@@ -386,18 +394,18 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
 
     from sumpy.visualization import FieldPlotter
 
-    eval_offset = np.array([5.5, 0.0])
+    eval_offset = np.array([5.5, 0.0, 0][:knl.dim])
 
     centers = (np.array(
             [
                 # box 0: particles, first mpole here
-                [0, 0],
+                [0, 0, 0][:knl.dim],
 
                 # box 1: second mpole here
-                np.array([-0.2, 0.1], np.float64),
+                np.array([-0.2, 0.1, 0][:knl.dim], np.float64),
 
                 # box 2: first local here
-                eval_offset + np.array([0.3, -0.2], np.float64),
+                eval_offset + np.array([0.3, -0.2, 0][:knl.dim], np.float64),
 
                 # box 3: second local and eval here
                 eval_offset