We consider the Dirichlet problem for a partial differential equation involving the Jacobian determinant in two dimensions of space. The problem consists in finding a vector-valued function such that the determinant of its gradient is given pointwise in a bounded domain, together with essential boundary conditions. This problem was initially considered in Dacorogna and Moser [Ann. Inst. H. Poincaré Anal. Non Linéaire, 7 (1990), pp. 1--26], and several theoretical generalizations have been derived since. In this work, we design a numerical algorithm for the approximation of the solution of such a problem for various kinds of boundary data. The proposed method relies on an augmented Lagrangian algorithm with biharmonic regularization, and low order mixed finite element approximations. An iterative method allows us to decouple the nonlinearity and the differential operators. Numerical experiments show the capabilities of the method for benchmarks and then for more demanding test problems.