Hyperbolic relaxation systems is an active field of research, with a largenumber of applications in physical modeling. Examples include modelsfor traffic flow, kinetic theory and fluid mechanics. This master’s thesis is a numerical and theoretical analysis of such systems, and consists of two main parts: The first is a new scheme for the stable numerical solution of hyperbolic relaxation systems using exponential integrators. First and second-order schemes of this type are derived and some desirable stability and accuracy properties are shown. The scheme is also used to solve a granular-gas model in order to demonstratethe practical use of the method. The second and largest part of this thesis is the analysis of the solutionsto 2 × 2 relaxation systems. In this work, the link between the the sub-characteristic condition and the stability of the solution of the relaxationsystem is discussed. In this context, the sub-characteristic condition andthe dissipativity of the Chapman–Enskog approximation are shown to beequivalent in both 1-D and 2-D. Also, the dispersive wave dynamics of hyperbolic relaxation systems isanalyzed in detail. For 2 × 2 systems, the wave-speeds of the individualFourier-components of the solution are shown to fulfill a transitional sub-characteristic condition. Moreover, the transition is monotonic in thevariable ξ = kε, where ε is the relaxation time of the system and k is thewave-number. A basic 2 × 2 model is used both as an example-model in the analyticaldiscussions, and as a model for numerical tests in order to demonstratethe implications of the analytical results.