The cross section for e^+e^- to eta J/psi between sqrt{s}=3.8 GeV/c^2 and 5.3 GeV/c^2 is measured via initial state radiation using 980 fb^{-1} of data on and around the Upsilon(nS)(n=1,2,3,4,5) resonances collected with the Belle detector at KEKB. Two resonant structures at the psi(4040) and psi(4160) are observed in the eta J/psi invariant mass distribution. Fitting the mass spectrum with the coherent sum of two Breit-Wigner functions, one obtains BR(psi(4040)toeta J/psi)cdotGamma_{ee}^{psi(4040)} = (4.8pm0.9pm1.4) eV and BR(psi(4160)toeta J/psi)cdotGamma_{ee}^{psi(4160)} = (4.0pm0.8pm1.4) eV for one solution and BR(psi(4040)toeta J/psi)cdotGamma_{ee}^{psi(4040)} = (11.2pm1.3pm1.9) eV and BR(psi(4160)toeta J/psi)cdotGamma_{ee}^{psi(4160)} = (13.8pm1.3pm2.0) eV for the other solution, where the first errors are statistical and the second are systematic. This is the first measurement of this hadronic transition mode of these two states, and the partial widths to eta J/psi are found to be about 1 MeV. There is no evidence for the Y(4260), Y(4360), psi(4415), or Y(4660) in the eta J/psi final state, and upper limits of their production rates in e^+e^- annihilation are determined.