Gravitational wave and electromagnetic observations can provide new insights into the nature of matter at supra-nuclear densities inside neutron stars. Improvements in electromagnetic and gravitational wave sensing instruments continue to enhance the accuracy with which they can measure the masses, radii, and tidal deformability of neutron stars. These better measurements place tighter constraints on the equation of state of cold matter above nuclear density. In this article, we discuss a complementary approach to get insights into the structure of neutron stars by providing a model prediction for non-linear fundamental eigenmodes (f-modes) and their decay over time, which are thought to be induced by time-dependent tides in neutron star binaries. Building on pioneering studies that relate the properties of f-modes to the structure of neutron stars, we systematically study this link in the non-perturbative regime using models that utilize numerical relativity. Using a suite of fully relativistic numerical relativity simulations of oscillating TOV stars, we establish blueprints for the numerical accuracy needed to accurately compute the frequency and damping times of f-mode oscillations, which we expect to be a good guide for the requirements in the binary case. We show that the resulting f-mode frequencies match established results from linear perturbation theory, but the damping times within numerical errors depart from linear predictions. This work lays the foundation for upcoming studies aimed at a comparison of theoretical models of f-mode signatures in gravitational waves, and their uncertainties with actual gravitational wave data, searching for neutron star binaries on highly eccentric orbits, and probing neutron star structure at high densities.