We revisit in one dimension the waterbag method to solve numerically Vlasov–Poisson equations. In this approach, the phase-space distribution function f (x, v) is initially sampled by an ensemble of patches, the waterbags, where f is assumed to be constant. As a consequence of Liouville theorem, it is only needed to follow the evolution of the border of these waterbags, which can be done by employing an orientated, self-adaptive polygon tracing isocontours of f. This method, which is entropy conserving in essence, is very accurate and can trace very well non-linear instabilities as illustrated by specific examples. As an application of the method, we generate an ensemble of single-waterbag simulations with decreasing thickness to perform a convergence study to the cold case. Our measurements show that the system relaxes to a steady state where the gravitational potential profile is a power law of slowly varying index β, with β close to 3/2 as found in the literature. However, detailed analysis of the properties of the gravitational potential shows that at the centre, β > 1.54. Moreover, our measurements are consistent with the value β = 8/5 = 1.6 that can be analytically derived by assuming that the average of the phase-space density per energy level obtained at crossing times is conserved during the mixing phase. These results are incompatible with the logarithmic slope of the projected density profile β − 2 ≃ −0.47 obtained recently by Schulz et al. using an N-body technique. This sheds again strong doubts on the capability of N-body techniques to converge to the correct steady state expected in the continuous limit.