Upgrades to existing gravitational-wave observatories have the potential to simultaneously constrain the nuclear equation of state and Hubble’s constant H₀ to per cent level with merging neutron star binaries. In practice, performing simultaneous inference of H₀ and the equation of state is limited computationally by the requirement to solve the equations of general-relativistic hydrostatic equilibrium millions of times. We develop a machine-learning model to solve the Tolman–Oppenheimer–Volkoff equations in less than a millisecond, and demonstrate its utility by performing direct inference of both equation of state and Hubble’s constant for synthetic neutron star merger signals with LIGO–Virgo–KAGRA operating at A+ sensitivities. We show that a population of 15 mergers observed with A+ allows for the radius of a 1.4 M_⊙ neutron star and H₀ to be constrained to km and at 90 per cent credible interval and 68 per cent credible interval, respectively. These constraints utilize only the gravitational-wave information to infer cosmological parameters; such numbers will be further improved with the addition of electromagnetic counterparts and/or galaxy catalogues.