We introduce a method for calculating the stationary state of a translation invariant array of weakly coupled cavities in the presence of dissipation and coherent as well as incoherent drives. Instead of computing the full density matrix our method directly calculates the correlation functions which are relevant for obtaining all local quantities of interest. It considers an expansion of the correlation functions and their equations of motion in powers of the photon tunneling rate between adjacent cavities, leading to an exact second order solution for any number of cavities. Our method provides a controllable approximation for weak tunneling rates applicable to the strongly correlated regime that is dominated by nonlinearities in the cavities and thus of high interest.