  We present the high-mass end of the galaxy stellar mass function using the
largest sample to date (5,352) of star-forming galaxies with $M_{\star} >
10^{11} M_{\odot}$ at cosmic noon, $1.5 < z < 3.5$. This sample is uniformly
selected across 17.2 deg$^2$ ($\sim$0.44 Gpc$^3$ comoving volume from $1.5 < z
< 3.5$), mitigating the effects of cosmic variance and encompassing a wide
range of environments. This area, a factor of 10 larger than previous studies,
provides robust statistics at the high-mass end. Using multi-wavelength data in
the Spitzer/HETDEX Exploratory Large Area (SHELA) footprint we find that the
SHELA footprint star-forming galaxy stellar mass function is steeply declining
at the high-mass end probing values as high as $\sim$$10^{-4}$ Mpc$^3$/dex and
as low as $\sim$5$\times$$10^{-8}$ Mpc$^3$/dex across a stellar mass range of
log($M_\star$/$M_\odot$) $\sim$ 11 - 12. We compare our empirical star-forming
galaxy stellar mass function at the high mass end to three types of numerical
models: hydrodynamical models from IllustrisTNG, abundance matching from the
UniverseMachine, and three different semi-analytic models (SAMs; SAG, SAGE,
GALACTICUS). At redshifts $1.5 < z < 3.5$ we find that results from
IllustrisTNG and abundance matching models agree within a factor of $\sim$2 to
10, however the three SAMs strongly underestimate (up to a factor of 1,000) the
number density of massive galaxies. We discuss the implications of these
results for our understanding of galaxy evolution.
