The paper presents an algorithm for obtaining the frequency responses of large-scale dynamical systems containing clearance type nonlinearities. The algorithm is based on coupling series of linear equations of motion, associated with piecewise linear functions. The accuracy of the solution does not significant depend on a magnitude of the integration step since the local linear equations of motion are solved analytically. The time step should be only sufficiently small for a reliable numerical determination of the switching point between two local linear equations of motion. It is a remarkable advantage with respect to other explicit integration methods. The proposed algorithm is validated in numerical example of the three-degree-of-freedom semi-definite system with two clearances under harmonic excitations. Good correlations with the results obtained by the describing function method, time finite element method and numerical integration confirm the algorithm.