Magnetohydrodynamics (MHD) is the physics branch that studies electrically conducting fluids under external magnetic fields. This thesis deals with the numerical approximation using stabilized finite element methods of two different formulations to model incompressible MHD, namely the resistive and inductionless MHD problems. Further, the linear systems of equations resulting from the application of these discrete formulations to simulate real cases are typically ill-conditioned and can have as many as 10̂6-10̂9 degrees of freedom. An efficient and scalable solver strategy is mandatory in these cases. On one hand, a new stabilized finite element formulation for the approximation of the resistive magnetohydrodynamics equations has been proposed. The novelty of this formulation with respect to existing ones is that it always converges to the physical solution, even when it is singular, which has been proved through a detailed stability and convergence analysis of the formulation. Moreover, it is inferred from the convergence analysis that a particular type of meshes with a macro-element structure is needed, which can be easily obtained after a straight modification of any original mesh. Finally, different operator splitting schemes have been proposed for solving the transient incompressible resistive MHD system that are unconditionally stable. Two levels of splitting have been considered. On the first level, the segregation of the Lagrange multipliers, the fluid pressure and the magneric pseudo-pressure, from the vectorial fields computation is achieved. On the second level, the fluid velocity and induction fields are also decoupled. This way, the fully coupled indefinite multiphysics system is transformed into smaller uncoupled one-physics problems. On the other hand, a stabilized formulation to solve the inductionless magnetohydrodynamic problem using the finite element method is presented. The inductionless MHD problem models the flow of an electrically charged fluid under the influence of an external magnetic field where the magnetic field induced in the fluid by the currents is negligible with respect to the external one. This system of partial differential equations is strongly coupled and highly nonlinear for real cases of interest. Therefore, solving the multiphysics linear systems of equations resulting from the discretization of these equations with finite element methods is a very challenging task which requires efficient and scalable preconditioners. A new family of recursive block LU preconditioners has been designed to improve the convergence of iterative solvers for this problem. These preconditioners are obtained after splitting the fully coupled matrix into one-physics problems for every variable (velocity, pressure, current density and electric potential) that can be optimally solved, e.g. using preconditioned domain decomposition algorithms. Furthermore, these ideas have been extended for developing recursive block LU preconditioners for the thermally coupled inductionless MHD problem.