diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index c060b8da9b47f2456580acad30aab5faa5fe7004..d43379a6a3e7e0fefa8e6fac2b93b36a7dd9070e 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -408,9 +408,9 @@ class NewLinearRecurrenceBasedDerivativeWrangler \
         if len(pde_mat) > 0:
             pde_mat = np.array(pde_mat, dtype=np.float64)
             n = nullspace(pde_mat, atol=tol)
-            k, idx, _ = interp_decomp(n.T, tol)
+            k, idx, proj = interp_decomp(n.T, tol)
+            s = np.hstack([np.eye(k), proj])[:, np.argsort(idx)]
             idx = idx[:k]
-            s = np.linalg.solve(n[idx, :].T, n.T).T
             stored_identifiers = [mis[i] for i in idx]
         else:
             s = np.eye(len(mis))