The problem of finding a $L^∞$-bounded two-dimensional vector field whose divergence is given in $L^2$ is discussed from the numerical viewpoint. A systematic way to find such a vector field is to introduce a non-smooth variational problem involving a $L^∞$-norm. To solve this problem from calculus of variations, we use a method relying on a well-chosen augmented Lagrangian functional and on a mixed finite element approximation. An Uzawa algorithm allows to decouple the differential operators from the nonlinearities introduced by the $L^∞$-norm, and leads to the solution of a sequence of Stokes-like systems and of an infinite family of local nonlinear problems. A simpler method, based on a $L^2$-regularization is also considered. Numerical experiments are performed, making use of appropriate numerical integration techniques when non-smooth data are considered; they allow to compare the merits of the two approaches discussed in this article and to show the ability of the related methods at capturing $L^∞$-bounded solutions.