In this paper, we consider the electromagnetically and thermally driven flow which is modeled by evolutionary magnetohydrodynamic equations and heat equation coupled through generalized Boussinesq approximation with temperature-dependent coefficients. Based on a third-order backward differential formula for temporal discretization, mixed finite element approximation for spatial discretization and extrapolated treatments in linearization for nonlinear terms, a linearized backward differentiation formula type scheme for the considered equations is proposed and analysed. Optimal $L^2$-error estimates for the proposed fully discretized scheme are obtained by the temporal-spatial error splitting technique. Numerical examples are presented to check the accuracy and efficiency of the scheme.