The dynamics of a one-dimensional quantum system coupled to a harmonic bath can be expressed through Feynman＇s path integral expression for the reduced density matrix. In this expression the influence of the environment is seen in correlations between positions of the system that are nonlocal in time. Makri and Makarov [J. Chem. Phys. 102, 4600 (1995)] showed that for many practical problems correlations over only a few time steps, Delta k(max) need to be taken into account, which led to an efficient iterative scheme. However, this algorithm scales as the size of the system to the power of 2(Delta k(max)+1), which restricts the size of the system that can be studied with this method. In this work we present an efficient algorithm which scales linearly with Delta k(max). In our method the reduced density matrix is written as a convolution of its past values with an integral equation kernel. The calculation of that kernel is based on a perturbative expansion of the discretized quasiadiabatic path integral expression for the reduced density matrix. The expansion ignores certain types of correlations.