diff --git a/devito/types/petsc.py b/devito/types/petsc.py
index c9307d2a09..f01eba7b65 100644
--- a/devito/types/petsc.py
+++ b/devito/types/petsc.py
@@ -5,6 +5,7 @@
 from cached_property import cached_property
 from devito.finite_differences import Differentiable
 from devito.types.basic import AbstractFunction
+from devito.finite_differences.tools import fd_weights_registry
 
 
 class DM(LocalObject):
@@ -72,13 +73,31 @@ class PETScArray(ArrayBasic, Differentiable):
 
     _data_alignment = False
 
+    # Default method for the finite difference approximation weights computation.
+    _default_fd = 'taylor'
+
     __rkwargs__ = (AbstractFunction.__rkwargs__ +
-                   ('dimensions', 'liveness'))
+                   ('dimensions', 'liveness', 'coefficients'))
+
+    def __init_finalize__(self, *args, **kwargs):
+
+        super().__init_finalize__(*args, **kwargs)
+
+        # Symbolic (finite difference) coefficients
+        self._coefficients = kwargs.get('coefficients', self._default_fd)
+        if self._coefficients not in fd_weights_registry:
+            raise ValueError("coefficients must be one of %s"
+                             " not %s" % (str(fd_weights_registry), self._coefficients))
 
     @classmethod
     def __dtype_setup__(cls, **kwargs):
         return kwargs.get('dtype', np.float32)
 
+    @property
+    def coefficients(self):
+        """Form of the coefficients of the function."""
+        return self._coefficients
+
     @cached_property
     def _C_ctype(self):
         petsc_type = dtype_to_petsctype(self.dtype)
diff --git a/tests/test_petsc.py b/tests/test_petsc.py
index ce3873b70f..955f361f7f 100644
--- a/tests/test_petsc.py
+++ b/tests/test_petsc.py
@@ -84,3 +84,7 @@ def test_petsc_subs():
 
     assert str(eqn_subs) == 'Eq(f1(x, y), Derivative(arr(x, y), (x, 2))' +  \
         ' + Derivative(arr(x, y), (y, 2)))'
+
+    assert str(eqn_subs.rhs.evaluate) == '-2.0*arr(x, y)/h_x**2' + \
+        ' + arr(x - h_x, y)/h_x**2 + arr(x + h_x, y)/h_x**2 - 2.0*arr(x, y)/h_y**2' + \
+        ' + arr(x, y - h_y)/h_y**2 + arr(x, y + h_y)/h_y**2'