Fluid systems with different densities are used in a wide range of applications in science and engineering, including geophysical and chemical flows. Numerical approximation of such problems can become cumbersome and computationally expensive, particularly when using classical monolithic schemes to approximate the Navier-Stokes system for incompressible flow. For this reason, we present a novel methodology for an unconditionally stable fractional step approach that treats the viscous term in an implicit-explicit (IMEX) manner, to linearize the Navier-Stokes system and reduce computational load. Additionally, the proposed method allows for the separation of velocity components, facilitating the solution of the momentum equations. We validate our proposed IMEX method through numerical experiments using finite element methods for spatial discretization. The tests range from manufactured solutions to complex two-phase Viscoplastic flows. We achieve the expected convergence order for each manufactured problem: at least O(∆t) for the fractional-step method, demonstrating both accuracy and stability. Additionally, we benchmark our method against a popular variable-viscosity test case, the Rayleigh-Taylor instability [1]. The evolution of the problem and the height of the rising bubble show good agreement with results reported by [1, 2]. Finally, we extend our method to the classical falling droplet benchmark [3], incorporating a Viscoplastic fluid using a regularized Bingham model. The results exhibit correct physical behavior without any numerical instabilities, highlighting the robustness of our method in handling complex physical phenomena without the need for additional stabilization techniques.