The close-packed AB$_2$ structures called Laves phases constitute the largest group of intermetallic compounds. In this paper we computationally investigated the pseudo-binary Laves phase system Y$_{1-x}$Gd$_x$(Fe$_{1-y}$Co$_y$)$_2$ spanning between the YFe$_2$, YCo$_2$, GdFe$_2$, and GdCo$_2$ vertices. While the vast majority of the Y$_{1-x}$Gd$_x$(Fe$_{1-y}$Co$_y$)$_2$ phase diagram is the ferrimagnetic phase, YCo$_2$ along with a narrow range of concentrations around it is the paramagnetic phase. We presented results obtained by Monte Carlo simulations of the Heisenberg model with parameters derived from first-principles calculations. For calculations, we used the Uppsala atomistic spin dynamics (UppASD) code together with the spin-polarized relativistic Korringa-Kohn-Rostoker (SPR-KKR) code. From first principles we calculated the magnetic moments and exchange integrals for the considered pseudo-binary system, together with spin-polarized densities of states for boundary compositions. Furthermore, we showed how the compensation point with the effective zero total moment depends on the concentration in the considered ferrimagnetic phases. However, the main result of our study was the determination of the Curie temperature dependence for the system Y$_{1-x}$Gd$_x$(Fe$_{1-y}$Co$_y$)$_2$. Except for the paramagnetic region around YCo$_2$, the predicted temperatures were in good qualitative and quantitative agreement with experimental results, which confirmed the ability of the method to predict magnetic transition temperatures for systems containing up to three different magnetic elements (Fe, Co, and Gd) simultaneously. For the Y(Fe$_{1-y}$Co$_y$)$_2$ and Gd(Fe$_{1-y}$Co$_y$)$_2$ systems our calculations matched the experimentally-confirmed Slater-Pauling-like behavior of T$_C$ dependence on the Co concentration.