diff --git a/sumpy/expansion/diff_op.py b/sumpy/expansion/diff_op.py
index db8e97e8a350054c315312837488c6c337d1b9e2..2a028859becf7e478811c602bb05ace3df5d1468 100644
--- a/sumpy/expansion/diff_op.py
+++ b/sumpy/expansion/diff_op.py
@@ -65,6 +65,13 @@ class LinearPDESystemOperator:
     def __hash__(self):
         return hash((self.dim, self.eqs))
 
+    @property
+    def order(self):
+        deg = 0
+        for eq in self.eqs:
+            deg = max(deg, max(sum(ident.mi) for ident in eq.keys()))
+        return deg
+
     def __mul__(self, param):
         eqs = []
         for eq in self.eqs:
diff --git a/test/test_misc.py b/test/test_misc.py
index 3a2055fa9082875aa9c153b691d57459358b0327..b0db855feabc11e43ad1e20590b3964bf8c95400 100644
--- a/test/test_misc.py
+++ b/test/test_misc.py
@@ -411,8 +411,11 @@ def test_as_scalar_pde_elasticity():
     ]
 
     pde = concat(*exprs)
+    assert pde.order == 1
     for i in range(5):
-        assert as_scalar_pde(pde, i) == laplacian(laplacian(diff_op[0]))
+        scalar_pde = as_scalar_pde(pde, i)
+        assert scalar_pde == laplacian(laplacian(diff_op[0]))
+        assert scalar_pde.order == 4
 
 
 def test_elasticity_new():