We present a novel approach for high-order accurate numerical differentiation on unstructured meshes of quadrilateral elements. To differentiate a given function, an auxiliary function with greater smoothness properties is defined which when differentiated provides the derivatives of the original function. The method generalises traditional finite difference methods to meshes of arbitrary topology in any number of dimensions for any order of derivative and accuracy. We confirm the accuracy of the numerical scheme using dual quadrilateral meshes and a refinement method based on subdivision surfaces. The scheme is applied to the solution of a range of partial differential equations, including both linear and nonlinear, and second and fourth order equations. For convection-dominated problems, we stabilise the discretisation using a generalised upwinding technique. Using relevant Navier-Stokes model problems, we compare our scheme to other high-order discretisations such as the discontinuous Galerkin method and demonstrate very high performance gains for equal accuracy.