The propagation of N photons in one dimensional waveguides coupled to M qubits is discussed, both in the strong and ultrastrong qubit–waveguide coupling. Special emphasis is placed on the characterisation of the nonlinear response and its linear limit for the scattered photons as a function of N, M, qubit inter distance and light–matter coupling. The quantum evolution is numerically solved via the matrix product states technique. The time evolutions for both the field and qubits are computed. The nonlinear character (as a function of N/M) depends on the computed observable. While perfect reflection is obtained for N/M ≅ 1, photon–photon correlations are still resolved for ratios N/M = non-zero. Inter-qubit distance enhances the nonlinear response. Moving to the ultrastrong coupling regime, we observe that inelastic processes are robust against the number of qubits and that the qubit–qubit interaction mediated by the photons is qualitatively modified. The theory developed in this work models experiments in circuit QED, photonic crystals and dielectric waveguides.