A high-fidelity finite volume scheme based on the BVD (boundary variation diminishing) concept is proposed in this study to solve the ideal magnetohydrodynamics (MHD) equations. A hybrid spatial reconstruction profile, consisting of a quadratic polynomial and a steepness-adjustable hyperbolic tangent function, is adopted to reproduce the accurate solutions of the complex magnetohydrodynamics flows. The BVD principle is used to find a optimal combination of these two types of spatial reconstructions by comparing the variations of the interface values interpolated in two adjacent cells, aiming to remove the non-physical oscillations around discontinuities by switching the quadratic polynomial to a step-shaped function. Additionally, a constrained transport (CT) method is applied in this study to assure the non-divergent solution of the magnetic field. The widely used numerical tests in one and two dimensional cases were checked in this study. The numerical results can retrieve the accuracy of a $3^{rd}$-order linear scheme in the convergence test for both advection and MHD equations and capture the strong shock waves in MHD flows without spurious oscillations. In comparison with the results of a $3^{rd}$-order WENO (weighted essentially non-oscillatory) scheme, the proposed scheme gains more accurate solutions not only for the strong discontinuities but also the smooth structures across scales. To our knowledge, this is the first attempt to build a high-fidelity model for ideal MHD equations by a BVD algorithm. Numerical results are competitive to those of existing advanced models and thus the BVD algorithm has promising potentials to build practical models for various MHD flows.