The applicability of nuclear ab initio calculations has rapidly extended over the past decades. However, starting research projects is still challenging due to the required numerical expertise in the generation of underlying nuclear interaction matrix elements and many-body calculations. To ease the first issue, in this paper we introduce the numerical code NuHamil to generate the nucleon-nucleon (NN) and three-nucleon (3N) matrix elements expressed in a spherical harmonic-oscillator basis, inputs of many-body calculations. The ground-state energies for the selected doubly closed shell nuclei are calculated with the no-core shell-model (NCSM) and in-medium similarity renormalization group (IMSRG). The code is written in modern Fortran, and OpenMP+MPI hybrid parallelization is available for the 3N matrix-element calculations.