  We introduce a general numerical method to compute dynamics and multi-time
correlations of chains of quantum systems, where each system may couple
strongly to a structured environment. The method combines the process tensor
formalism for general (possibly non-Markovian) open quantum systems with time
evolving block decimation (TEBD) for 1D chains. It systematically reduces the
numerical complexity originating from system-environment correlations before
integrating them into the full many-body problem, making a wide range of
applications numerically feasible. We illustrate the power of this method by
studying two examples. First, we study the thermalization of individual spins
of a short XYZ Heisenberg chain with strongly coupled thermal leads. Our
results confirm the complete thermalization of the chain when coupled to a
single bath, and reveal distinct effective temperatures in low, mid, and high
frequency regimes when the chain is placed between a hot and a cold bath.
Second, we study the dynamics of diffusion in an longer XY chain, when each
site couples to its own bath.
