Efficient computer codes for the explicitly correlated coupled-cluster (CC-R12 or F12) methods with up to triple (CCSDT-R12) and quadruple excitations (CCSDTQ-R12), which take account of the spin, Abelian point-group, and index-permutation symmetries and are based on complete diagrammatic equations, have been implemented with the aid of the computerized symbolic algebra SMITH. Together with the explicitly correlated coupled-cluster singles and doubles (CCSD-R12) method reported earlier [T. Shiozaki et al., J. Chem. Phys. 129, 071101 (2008)], they form a hierarchy of systematic approximations (CCSD-R12<CCSDT-R12<CCSDTQ-R12) that converge very rapidly toward the exact solutions of the polyatomic Schrodinger equations with respect to both the highest excitation rank and basis-set size. Using the Slater-type function exp(-gamma r(12)) as a correlation function, a CC-R12 method can provide the aug-cc-pV5Z-quality results of the conventional CC method of the same excitation rank using only the aug-cc-pVTZ basis set. Combining these CC-R12 methods with the grid-based, numerical Hartree-Fock equation solver [T. Shiozaki and S. Hirata, Phys. Rev. A 76, 040503(R) (2007)], the solutions (eigenvalues) of the Schrodinger equations of neon, boron hydride, hydrogen fluoride, and water at their equilibrium geometries have been obtained as -128.9377+/-0.0004, -25.2892+/-0.0002, -100.459+/-0.001, and -76.437+/-0.003 E(h), respectively, without resorting to complete-basis-set extrapolations. These absolute total energies or the corresponding correlation energies agree within the quoted uncertainty with the accurate, nonrelativistic, Born-Oppenheimer values derived experimentally and/or computationally.