  Magnetic order in two-dimensional (2D) materials is intimately coupled to
magnetic anisotropy (MA) since the Mermin-Wagner theorem implies that
rotational symmetry cannot be spontaneously broken at finite temperatures in
2D. Large MA thus comprises a key ingredient in the search for magnetic 2D
materials that retains the magnetic order above room temperature. Magnetic
interactions are typically modeled in terms of Heisenberg models and the
temperature dependence on magnetic properties can be obtained with the Random
Phase Approximation (RPA), which treats magnon interactions at the mean-field
level. In the present work we show that large MA gives rise to strong
magnon-magnon interactions that leads to a drastic failure of the RPA. We then
demonstrate that classical Monte Carlo (MC) simulations correctly describe the
critical temperatures in the large MA limit and agree with RPA when the MA
becomes small. A fit of the MC results leads to a simple expression for the
critical temperatures as a function of MA and exchange coupling constants,
which significantly simplifies the theoretical search for new 2D magnetic
materials with high critical temperatures. The expression is tested on a
monolayer of CrI$_3$, which were recently observed to exhibit ferromagnetic
order below 45 K and we find excellent agreement with the experimental value.
