In the field of Computational Fluid Dynamics (CFD), solving conjugate heat transfer problems involving natural convection remains a computationally demanding endeavor. This study presents the NEMESYS algorithm, specifically optimized for Graphics Processing Units (GPUs), with an application focused on natural convection in a square cavity that incorporates a centrally situated solid block. This complex configuration necessitates the concurrent resolution of the Navier-Stokes equations governing fluid flow and the energy equation governing heat conduction in both fluid and solid phases. Through the efficient utilization of GPU-based parallel processing, the NEMESYS algorithm managed to markedly reduce the computational burden. A quantifiable 99.5% time reduction was recorded when compared to equivalent Central Processing Unit (CPU)-based simulations, thereby manifesting a significant leap in computational efficiency. To authenticate the algorithm's credibility, an exhaustive validation process was undertaken. The simulation results were cross-verified against established benchmarks from academic literature as well as outputs from widely used commercial CFD software packages. This validation revealed strong agreement in critical parameters such as fluid velocity and temperature distributions within the fluid cavity, as well as heat conduction characteristics within the solid block. In summary, the NEMESYS algorithm emerges as a reliable, efficient computational framework for tackling the intricacies of natural convection problems involving conjugate heat transfer and holds potential for broader adaptability in high-fidelity CFD simulations.