Feedback from massive stars is believed to play a critical role in shaping the galaxy mass function, the structure of the interstellar medium (ISM) and the low efficiency of star formation, but the exact form of the feedback is uncertain. In this paper, the first in a series, we present and test a novel numerical implementation of stellar feedback resulting from momentum imparted to the ISM by radiation, supernovae and stellar winds. We employ a realistic cooling function, and find that a large fraction of the gas cools to ≲100 K, so that the ISM becomes highly inhomogeneous. Despite this, our simulated galaxies reach an approximate steady state, in which gas gravitationally collapses to form giant 'molecular' clouds (GMCs), dense clumps and stars; subsequently, stellar feedback disperses the GMCs, repopulating the diffuse ISM. This collapse and dispersal cycle is seen in models of Small Magellanic Cloud (SMC)-like dwarfs, the Milky Way and z∼ 2 clumpy disc analogues. The simulated global star formation efficiencies are consistent with the observed Kennicutt-Schmidt relation. Moreover, the star formation rates are nearly independent of the numerically imposed high-density star formation efficiency, density threshold and density scaling. This is a consequence of the fact that, in our simulations, star formation is regulated by stellar feedback limiting the amount of very dense gas available for forming stars. In contrast, in simulations without stellar feedback, i.e. under the action of only gravity and gravitationally induced turbulence, the ISM experiences runaway collapse to very high densities. In these simulations without feedback, the global star formation rates exceed observed galactic star formation rates by 1-2 orders of magnitude, demonstrating that stellar feedback is crucial to the regulation of star formation in galaxies.