We study the behavior of stationary non-equilibrium two-body correlation functions for Diffusive Systems with equilibrium reference states (DSe). A DSe is described at the mesoscopic level by $M$ locally conserved continuum fields that evolve through coupled Langevin equations with white noises. The dynamic is designed such that the system may reach equilibrium states for a set of boundary conditions. In this form, just by changing the equilibrium boundary conditions, we make the system driven to a non-equilibrium stationary state. We decompose the correlations in a known local equilibrium part and another one that contains the non-equilibrium behavior and that we call {it correlations excess} $bar C(x,z)$. We formally derive the differential equations for $bar C$. We define a perturbative expansion around the equilibrium state to solve them order by order. We show that the $bar C$s first-order expansion, $bar C^{(1)}$, is always zero for the unique field case, $M=1$. Moreover $bar C^{(1)}$ is always long-range or zero when $M>1$. Surprisingly we show that their associated fluctuations, the space integrals of $bar C^{(1)}$, are always zero. Therefore, the fluctuations are dominated by the local equilibrium behavior up to second order in the perturbative expansion around the equilibrium. We derive the behaviors of $bar C^{(1)}$ in real space for dimensions $d=1$ and $2$ explicitly, and we apply the analysis to a generic $M=2$ case and, in particular, to a hydrodynamic model where we explicitly compute the two first perturbative orders, $bar C^{(1),(2)}$, and its associated fluctuations.