An iterative algorithm is presented for evaluating the path integral expression for the reduced density matrix of a quantum system interacting with an anharmonic dissipative bath whose influence functional is obtained via numerical methods. The method allows calculation of the reduced density matrix over very long time periods.