Relativistic jets are intrinsic phenomena of active galactic nuclei (AGN) and quasars. They have been observed to also emanate from systems containing compact objects, such as white dwarfs, neutron stars and black hole candidates. The corresponding Lorentz factors, $Gamma$, were found to correlate with the compactness of the central objects. In the case of quasars and AGNs, plasmas with $Gamma$-factors larger than $8$ were detected. However, numerically consistent modelling of propagating shock-fronts with $Gamma geq 4$ is a difficult issue, as the non-linearities underlying the transport operators increase dramatically with $Gamma$, thereby giving rise to a numerical stagnation of the time-advancement procedure or alternatively they may diverge completely. In this paper, we present a unified numerical solver for modelling the propagation of one-dimensional shock fronts with high Lorentz factors. The numerical scheme is based on the finite-volume formulation with adaptive mesh refinement (AMR) and domain decomposition for parallel computation. It unifies both time-explicit and time-implicit numerical schemes within the framework of the pre-conditioned defect-correction iteration solution procedure. We find that time-implicit solution procedures are remarkably superior over their time-explicit counterparts in the very high $Gamma$-regime and therefore most suitable for consistent modelling of relativistic outflows in AGNs and micro-quasars.