Abstract The present work deals with the nonlinear multiphysics advection-diffusion problem where non-stationary Navier-Stokes equations for incompressible flows are coupled with heat transfer equations.Although a large variety of numerical methods is widespread in literature, a sufficiently comprehensive investigation, to the best of authors' knowledge, has not been provided yet when aiming for high performances on both accuracy and efficiency, either by solving these problems jointly or separately. To this end, we propose a modified Newton-Raphson (mNR) scheme, based on a fully nonlinear formulation and employing standard Lagrange Finite Elements, implemented in an in-house C++ code using the open-source deal.II library. The proposed method is first validated against consolidated results available in the literature, testing fluid motion equations first alone, and then coupled with thermal convections. Then, its performance is highlighted by comparison with enhanced schemes, ad-hoc implemented, approaching the coupled problems through projection methods. Through a wide numerical campaign, we finally prove that the mNR scheme is able to achieve higher level of accuracy and efficiency with respect to the other considered methods. The C++ code implementing the proposed mNR scheme together with those employed for comparisons is freely available at https://doi.org/10.5281/zenodo.8192391.