Plasma rotation proves to have important effects on the improvement of the confinement and on the entrance into the H-mode of improved confinement in tokamaks. To consider this issue, we propose to determine numerically the steady states of the visco-resistive MHD equations including the non-linear (v.∇)v term. We take into account the external electric field in the toroidal direction used to create the toroidal current required in a tokamak to create the poloidal component of the magnetic field. Numerical results in JET and ITER geometries obtained using the finite element programming language FreeFem++ are presented. We consider realistic values of the resistivity η and vary the viscosity, ν, of which the realistic order of magnitude is poorly known. The axisymmetric steady-state visco-resistive MHD equations with symmetric boundary conditions are solved using a continuation method on ν and a Newton-Raphson scheme to handle the nonlinearity.