In biological systems, liquid–liquid phase separation (LLPS) often occurs in confined and heterogeneous environments, such as cytoplasmic regions adjacent to lipid membranes. Capturing the coupled dynamics between bulk phase separation and membrane-associated processes therefore requires models and numerical methods that can faithfully represent both domains and their interactions. To address this challenge, we developed a mass-conserving, linear, and bound-preserving numerical scheme for a newly proposed bulk–membrane coupled phase-field model. Numerical experiments confirm the second-order temporal accuracy of the proposed scheme and demonstrate its ability to capture key biologically relevant phenomena, including droplet spreading kinetics on membranes with surface binding and phase separation driven by strong repulsive interactions among multiple components. Overall, this work provides a robust computational framework for investigating emergent behaviors in bulk–membrane coupled systems under biologically relevant conditions. The proposed model and numerical scheme offer new quantitative insights into the spatial organization of membrane-associated biomolecular condensates and advance our understanding of membrane-mediated phase behavior in living cells.