Complete information on the equilibrium behaviour and dynamics of a quantum field theory (QFT) is provided by multipoint correlation functions. However, their theoretical calculation is a challenging problem, even for exactly solvable models. This has recently become an experimentally relevant problem, due to progress in cold-atom experiments simulating QFT models and directly measuring higher order correlations. Here we compute correlation functions of the quantum sine-Gordon model, a prototype integrable model of central interest from both theoretical and experimental points of view. Building upon the so-called Truncated Conformal Space Approach, we numerically construct higher order correlations in a system of finite size in various physical states of experimental relevance, both in and out of equilibrium. We measure deviations from Gaussianity due to the presence of interaction and analyse their dependence on temperature, explaining the experimentally observed crossover between Gaussian and non-Gaussian regimes. We find that correlations of excited states are markedly different from the thermal case, which can be explained by the integrability of the system. We also study dynamics after a quench, observing the effects of the interaction on the time evolution of correlation functions, their spatial dependence, and their non-Gaussianity as measured by the kurtosis.