Anisotropic thermal conduction plays an important role in determining the structure of the hot plasma in the solar corona. When hot plasma appears, the conductivity rises with temperature and becomes highly nonlinear. Explicit solvers for parabolic problems often lead to much smaller time-steps limited by a Courant–Friedrichs–Lewy (CFL) condition in comparison with hyperbolic Magnetohydrodynamics (MHD) equations. In this work, we present a pseudo-linear, directionally-split, semi-implicit method allowing for large time-steps as well as the optimized parallelization algorithm, integrated with the MHD solver. Our scheme can perfectly preserve the monotonicity and the geometry of shocks and discontinuities in complex MHD problems.Two sets of numerical tests show that an increase in time step of $\sim 600$ can be easily achieved with an acceptable error by our scheme compared to explicit methods, and the use of large time-steps can still follow fast dynamic processes reliably. In addition, the extendibility studies have proven that the associated parallel efficiency is comparably high. This method is also useful for any kind of time-dependent conductivity problems for the solar applications in the future.