Solving partial differential equations (PDEs) in domains with complex geometries is essential for many applications but presents substantial computational challenges, especially when it comes to accurately imposing boundary conditions. This talk presents high-accuracy discretization techniques for boundary conditions, aimed at enhancing precision within the framework of unfitted boundary methods. We propose advanced high-order discretization schemes based on the ghost point method. This approach extends the computational grid beyond the physical domain by introducing ghost points, whose values are carefully computed to enforce boundary conditions with high fidelity. Unlike traditional methods that often assign ghost point values through simple extrapolation (independently from nearby points) our method employs a coupled approach. Here, ghost point values are solved alongside neighboring internal and other ghost points, forming an augmented linear system that improves both accuracy and stability. To solve the resulting augmented systems efficiently, we design a tailored multigrid solver optimized for curved boundaries. We validate the effectiveness of our approach through numerical experiments on elliptic PDEs and demonstrate its applicability to incompressible fluid simulations, including dynamic scenarios like oscillating bubbles. Although our method is implemented within a finite difference framework, its core principles extend naturally to finite volume and finite element formulations.