AbstractTransient responses of linear systems governed by partial dffferential equations (withtime and one spatial coordinate as independent variables), having time-dependent coefficients and non-homogeneous initial conditions,are simulated numerically by a hybrid use of Crump's scheme for the inversion of Laplace transforms and some traditional numerical methods. To handle the time dependency of coefficients, the“stairway”approximation whereby a time-curve is replaced by a series of small steps is incorporated. For each of these steps, the technique of Laplace transforms is applied; with an“initial”condition supplied by the solution at the end of the preceding step. These non-homogeneous“initial”conditions are fitted with high-degree polynomials so that the particular integrals to the ordinary differential equations in the trànsformed domain can be written down explicitly as polynomials. The solutions in the transformed domain can be, then, inverted using a numerical inversion scheme of the new generation (based on a Fourier series expansion). The hybrid method is applied,as an example, to the engineering system of packed-bed energy storage and regeneration where the mass-flow rateof fluid stream is time-dependent, and may change itssign (direction of flow). Numerical results show the method to be fast, robust, easy to implement, and with acceptabledegree of accuracy.