  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.
