<p style='text-indent:20px;'>This paper is devoted to the theoretical and numerical analysis of the non-stationary Rayleigh-Bénard-Marangoni (RBM) system. We analyze the existence of global weak solutions for the non-stationary RBM system in polyhedral domains of <inline-formula><tex-math id="M1">\begin{document}$ \mathbb{R}^3 $\end{document}</tex-math></inline-formula> and the convergence, in the norm of <inline-formula><tex-math id="M2">\begin{document}$ L^{2}(\Omega), $\end{document}</tex-math></inline-formula> to the corresponding stationary solution. Additionally, we develop a numerical scheme for approximating the weak solutions of the non-stationary RBM system, based on a mixed approximation: finite element approximation in space and finite differences in time. After proving the unconditional well-posedness of the numerical scheme, we analyze some error estimates and establish a convergence analysis. Finally, we present some numerical simulations to validate the behavior of our scheme.