There exist many methods to calculate forced response in mechanical systems. Some methods are slow and the errors introduced are unknown. The paper presents a method that uses digital filters and modal superposition. It is shown how aliasing can be avoided as well as phase errors. The parameters describing the mechanical system are residues and poles, taken from FEA models, from lumped MCK systems, from analytic solutions or from experimental modal analysis. Modal damping may be used. The error in the calculation is derived and is shown to be only a function of the sampling frequency used. When the method is applied to linear mechanical systems in MATLAB it is very fast. The method is extended to incorporate nonlinear components. The nonlinear components could be simple, like hardening or stiffening springs, but may also contain memory, like dampers with hysteresis. The simulations are used to generate test data for development and evaluation of methods for identification of non-linear systems.