  We investigate the impact of using high-order numerical methods to study the
merger of magnetised neutron stars with finite-temperature microphysics and
neutrino cooling in full general relativity. By implementing a fourth-order
accurate conservative finite-difference scheme we model the inspiral together
with the early post-merger and highlight the differences to traditional
second-order approaches at the various stages of the simulation. We find that
even for finite-temperature equations of state, convergence orders higher than
second order can be achieved in the inspiral and post-merger for the
gravitational-wave phase. We further demonstrate that the second-order scheme
overestimates the amount of proton-rich shock-heated ejecta, which can have an
impact on the modelling of the dynamical part of the kilonova emission.
Finally, we show that already at low resolution the growth rate of the magnetic
energy is consistently resolved by using a fourth-order scheme.
