In this paper, we consider the nonstationary magnetohydrodynamic coupled heat equation through the well-known Boussinesq approximation. The second order backward difference formula is used for time derivative terms, and the mixed finite method is used for spatial discretization, we employ the Taylor-Hood elements to approximate heat and Navier-Stokes equations, N$\mathrm{\acute{e}}$d$\mathrm{\acute{e}}$lec edge elements are used to approximate the magnetic induction. The divergence free conditions are weakly satisfied at the discrete level. Due to the use of N$\mathrm{\acute{e}}$d$\mathrm{\acute{e}}$lec edge element, the proposed method is particularly suitable for problems defined on non-smooth and multi-connected domains. Moreover, the numerical scheme is energy conserving. Under the weak regularity hypothesis of the exact solution, we present error estimate for velocity, magnetic variable and temperature. Finally, the convergence analysis is verified by some experiments, and the magnetic fluid phenomenon is simulated by driven cavity flow.