  We show how to exploit algebraic relations of operators (or matrices) which
constitute the non-equilibrium matrix product steady state of a boundary driven
quantum spin chain to find partial differential equations determining all the
$m$-point correlation functions in the continuum limit. These partial
differential equations, the order of which is determined by scaling of the
non-equilibrium partition function, are readily solved if we also know the
boundary conditions. In this way we can avoid resorting to representation
theory of the matrix product algebra. We apply our methods to study the
distributions, or moments, of the magnetization and the spin current
observables in boundary driven open XXZ spin chains, as well as some connected
correlation functions and find that the transverse connected correlation
functions are thermodynamically non-decaying and long-range at the isotropic
point $\Delta=1$.
