The solution of the elastohydrodynamic lubrication (EHL) problem involves the simultaneous resolution of the hydrodynamic (Reynolds equation) and elastic problems (elastic deformation of the contacting surfaces). Up to now, most of the numerical works dealing with the modeling of the isothermal EHL problem were based on a weak coupling resolution of the Reynolds and elasticity equations (semi-system approach). The latter were solved separately using iterative schemes and a finite difference discretization. Very few authors attempted to solve the problem in a fully coupled way, thus solving both equations simultaneously (full-system approach). These attempts suffered from a major drawback which is the almost full Jacobian matrix of the nonlinear system of equations. This work presents a new approach for solving the fully coupled isothermal elastohydrodynamic problem using a finite element discretization of the corresponding equations. The use of the finite element method allows the use of variable unstructured meshing and different types of elements within the same model which leads to a reduced size of the problem. The nonlinear system of equations is solved using a Newton procedure which provides faster convergence rates. Suitable stabilization techniques are used to extend the solution to the case of highly loaded contacts. The complexity is the same as for classical algorithms, but an improved convergence rate, a reduced size of the problem and a sparse Jacobian matrix are obtained. Thus, the computational effort, time and memory usage are considerably reduced.