We present a new method for the determination of highly precise continuous relaxation spectra. The method is based on the use of piecewise cubic Hermite splines, which are fairly easy to tabulate by using their knots. The Hermite splines method allows a continuous description of the spectrum by a series of polynomial functions. The numerical instabilities of the spectrum calculation are minimized by limiting the slope of the spectrum to physically meaningful values. The reproducibility of the spectrum calculation is within an error margin of about +/- 10% in the physically relevant relaxation time range. This method is able to retrieve the spectrum based on data calculated from a benchmark with high accuracy and precision.