  We constrain the space density and properties of massive galaxy candidates at
redshifts of z>=3.5 in the GOODS-N field. By selecting sources in the
Spitzer+IRAC bands, a highly stellar mass-complete sample is
assembled,including massive galaxies which are very faint in the
optical/near-IR bands that would be missed by samples selected at shorter
wavelengths. The z>=3.5 sample was selected down to 23 mag at 4.5 micron using
photometric redshifts that have been obtained by fitting the galaxies SEDs at
optical, near-IR and IRAC bands. We also require that the brightest band in
which candidates are detected is the IRAC 8 micron band in order to ensure that
the near-IR 1.6 micron (rest-frame) peak is falling in or beyond this band. We
found 53 z>=3.5 candidates, with masses in the range of M~10^10-10^11M_sun. At
least ~81% of these galaxies are missed by traditional Lyman Break selection
methods based on UV light. Spitzer+MIPS emission is detected for 60% of the
sample of z>=3.5 galaxy candidates. Although in some cases this might suggest a
residual contamination from lower redshift star-forming galaxies or AGN, 37% of
these objects are also detected in the sub-mm/mm bands in recent SCUBA,AzTEC
and MAMBO surveys, and have properties fully consistent with vigorous starburst
galaxies at z>=3.5. The comoving number density of galaxies with stellar masses
>= 5x10^10M_sun(a reasonable stellar mass completeness limit for our sample) is
2.6x10^-5Mpc^-3 (using the volume within 3.5<z<5), and the corresponding
stellar mass density is ~2.9x10^6M_sunMpc^-3, or~3% of the local density above
the same stellar mass limit.For the sub-sample of MIPS-undetected galaxies,we
find a number density of ~0.97x10^-5Mpc^-3 and a stellar mass density of
~1.15x10^6M_sun Mpc^-3.[abridged]
