The Navier-Stokes-Cahn-Hilliard (NSCH) diffuse interface model is a phase-field model that describes binary-fluid systems. This model is inherently capable of dealing with topological changes due to its phase-field nature and therefore particularly suitable for complex multi-phase flow. Recently, an advancement has been made to include complex domains using an immersed method (Stoter, Sluijs, Brummelen, Verhoosel, & Demont, 2023). We numerically approximate the NSCH model by means of an isogeometric analysis method with adaptive spatial resolution based on a residual-based error estimate. In addition, we use ε-continuation for robust spatial-refinements and to decouple the interface movement per time step from the interface width (Van Brummelen, Demont, & van Zwieten, 2020). In this model, the interface length scale is minuscule compared to the system length scale while encompassing the most complex dynamics. Therefore, this relatively small part of the domain carries a large part of the computational cost. To reduce the computational costs, we use high-order splines with a high continuity to capture the interface phenomena efficiently while limiting the number of degrees of freedom. To reduce the number of timesteps, we introduce a third-order generalized-α method. This two-step higher-order time discretization method does not impose the need for increased spatial resolution. Therefore, in this work, we present higher-order two-step time discretization methods for the NSCH model based on the generalized-α method as presented in (Behnoudfar, Deng, & Calo, 2020).