Multi-scale problems are omnipresent in environmental and industrial processes, posing a challenge to classical numerical solvers given that the propagation speeds of information span several orders of magnitude. In hyperbolic systems, the absolute fastest wave speed remains finite, but in explicitly integrated numerical schemes, it determines the time step that ensures stability. This means that it may lead to vanishing time increments in the presence of fast processes caused for instance by high pressure or strong magnetic fields. Furthermore, it is known that upwind schemes introduce spurious numerical diffusion into the approximate solution: a problem that can only be partially mitigated by mesh refinement. Therefore, a common approach, also applied here, is to use implicit or semi-implicit time integrators combined with centered differences for spatial derivatives in implicitly treated systems, in order to obtain scale-independent artificial dissipation and stability under large time steps. As the evolution of the modeled variables is described by a nonlinear flux function, treating it fully or partially implicitly involves solving nonlinear systems. Depending on the problem, these systems can be large, coupled, ill-posed systems, for which solvers such as the Newton method may converge very slowly or not at all. To avoid nonlinear implicit systems, we apply Jin-Xin relaxation to the implicitly treated stiff flux terms. This leads to a linear flux structure in the obtained relaxation model. Therein, the nonlinearity of the stiff flux is transferred to an algebraic relaxation source. By splitting away the relaxation source, for each variable, a decoupled wave-type equation can be written, yielding a predicted solution at temporarily frozen wave speeds. Projecting onto the relaxation equilibrium manifold, taking the relaxation source into account, we get a prediction-correction scheme for the original problem. Work with Puppo and Iollo