  We devise an iterative scheme for numerically calculating dynamical two-point
correlation functions in integrable many-body systems, in the Eulerian scaling
limit. Expressions for these were originally derived in Ref. [1] by combining
the fluctuation-dissipation principle with generalized hydrodynamics.
Crucially, the scheme is able to address non-stationary, inhomogeneous
situations, when motion occurs at the Euler-scale of hydrodynamics. In such
situations, in interacting systems, the simple correlations due to fluid modes
propagating with the flow receive subtle corrections, which we test. Using our
scheme, we study the spreading of correlations in several integrable models
from inhomogeneous initial states. For the classical hard-rod model we compare
our results with Monte-Carlo simulations and observe excellent agreement at
long time-scales, thus providing the first demonstration of validity for the
expressions derived in Ref. [1]. We also observe the onset of the Euler-scale
limit for the dynamical correlations.
