For some loading applications, the normal modes approach to the viscoelastic relaxation of a spherical earth requires the use of spherical harmonics up to a high degree. Examples include postseismic deformation (internal loading) and sea level variations due to glacial isostatic adjustment (external loading). In the case of postseismic modelling, the convergence of the solution, given as a spherical harmonic expansion series, is directly dependent on loading depth and requires several thousands of terms for shallow earthquake sources. The particular structure of the analytical fundamental solutions used in normal mode techniques usually does not allow a straightforward calculation, since numerical problems can readily occur due to the stiffness of the matrices used in the propagation routines. Here we show a way of removing this stiffness problem by approximating the fundamental matrix solutions, followed by a rescaling procedure, in this way we can virtually go up to whatever harmonic degree is required.