A numerical formulation to solve the MHD problem with thermal coupling is presented in full detail. The distinctive feature of the method is the design of the stabilization terms, which serve several purposes. First, convective dominated flows in the Navier-Stokes and the heat equation can be dealt with. Second, there is no restriction in the choice of the interpolation spaces of all the variables and, finally, flows highly coupled with the magnetic field can be accounted for. Different aspects related to the design of the final fully discrete and linearized algorithm are also discussed.