We address the numerical solution of 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 point-wise in a bounded domain, together with essential boundary conditions. The proposed numerical algorithm relies on an augmented Lagrangian algorithm with biharmonic regularization, and low order mixed finite element approximations. An iterative method allows to decouple the local nonlinearities and the global variational problem that involves a biharmonic operator. Numerical experiments validate the proposed method.