Binary neutron star mergers are associated with a variety of observable phenomena in the gravitational and electromagnetic spectra and are of great importance in a number of different physical subjects, e.g. high energy and gravitational physics. In this thesis, we are investigating binary neutron star systems in the last milliseconds before and after their merger. In such systems gravity is strong and has to be described by Einstein’s full Theory of General Relativity. Because of the complexity of the governing equations of general relativity and relativistic hydrodynamics no analytical solutions exist. Thus, the usage of numerical methods is inevitable. Throughout the thesis we consider different configurations by varying the spin, the equation of state, and the mass-ratio. In particular, we present the first consistent, constraint solved simulations of spinning binary neutron stars in full general relativity and the highest mass ratios simulated to date. Additionally, new numerical methods were implemented in the existing BAM code, most notably a refluxing algorithm which ensures mass conservation across mesh refinement boundaries. This algorithm allowed us to perform the most accurate simulations of the gravitational collapse of a rotating neutron star. In addition to pure numerical waveform modeling, we used high-resolution simulations to validate an improved tidal effective-one-body model and show that the new formalism can predict the waveform accurately up to the moment of merger. Furthermore, the effective-one-body model predicts quasi-universal relations, which we found also in full general relativistic simulations for the inspiral and even in the postmerger phase.