There have been developed a number of algorithms for the determination of relaxation spectrum from viscoelastic data. When viscoelastic data are given by stress relaxation test, the relation between relaxation spectrum and relaxation modulus can be transformed to that of Laplace transform. Hence, calculation of relaxation spectrum from relaxation modulus is a problem of inverting Laplace transform. Among various mathematical methods for inverse Laplace transform, the Post–Widder formula has been chosen by a number of researchers. However, they did not solve the problem by use of advance numerical technic but used low-order approximations. We suggest a new numerical algorithm which can calculate higher order solutions. In principle, the order of the Post–Widder formula is not limited in our algorithm.