Binary neutron star coalescence lead to extremely violent merger events in the universe. Such events give rise to a variety of observable phenomena in both the gravitational and the electromagnetic channels. In this thesis we study such extreme events in the last milliseconds around the merger. Gravity and its coupling to other physics for such systems is described by Einsteins theory of general relativity. Since these are highly nonlinear coupled set of partial differential equations, no analytical solutions exist for generic and dynamical systems such as binary neutron stars in the strong-field regime. Therefore, the usage of numerical methods is inevitable. In this thesis we consider configurations of binary neutron stars with varying equation of state, spin, eccentricity and spin orientation. In particular we present the first numerical relativity simulations of highly eccentric binary neutron stars in full (3+1)D using consistent initial data, i.e., setups which are in agreement with the Einstein equations and with the equations of general relativistic hydrodynamics in equilibrium. Additionally, we also study precessing binary neutron star systems using full (3+1)D numerical relativity simulations with consistent initial data. Further, we investigate a new numerical approach for solving the equations of general relativistic hydrodynamics using the discontinuous Galerkin method. It combines the key aspects of finite volume and traditional finite element methods, for e.g., element locality and hp adaptivity besides the benefit of spectral convergence for smooth solutions. Discontinuous Galerkin methods ease the treatment of complex grid geometries and, allow simple and efficient parallelization.