We infer the gravitational potential of the Galactic disk by analysing the phase-space densities of 120 stellar samples in 40 spatially separate sub-regions of the solar neighbourhood, using Gaias second data release (DR2), in order to quantify spatially dependent systematic effects that bias this type of measurement. The gravitational potential was inferred under the assumption of a steady state in the framework of a Bayesian hierarchical model. We performed a joint fit of our stellar tracers three-dimensional velocity distribution, while fully accounting for the astrometric uncertainties of all stars. The inferred gravitational potential is compared, post-inference, to a model for the baryonic matter and halo dark matter components. We see an unexpected but clear trend for all 40 spatially separate sub-regions: Compared to the potential derived from the baryonic model, the inferred gravitational potential is significantly steeper close to the Galactic mid-plane (<60 pc), but flattens such that the two agree well at greater distances (~400 pc). The inferred potential implies a total matter density distribution that is highly concentrated to the Galactic mid-plane and decays quickly with height. Apart from this, there are discrepancies between stellar samples, implying spatially dependent systematic effects which are, at least in part, explained by substructures in the phase-space distributions. In terms of the inferred matter density distribution, the very low matter density that is inferred at greater heights is inconsistent with the observed scale height and matter distribution of the stellar disk, which cannot be explained by a misunderstood density of cold gas or other hidden mass. Our interpretation is that these results must be biased by a time-varying phase-space structure, possibly a breathing mode, that is large enough to affect all stellar samples in the same manner.