  The formalism of matrix product states is used to perform a numerical study
of 1+1 dimensional QED -- also known as the (massive) Schwinger model -- in the
presence of an external static `quark' and `antiquark'. We obtain a detailed
picture of the transition from the confining state at short interquark
distances to the broken-string `hadronized' state at large distances and this
for a wide range of couplings, recovering the predicted behavior both in the
weak and strong coupling limit of the continuum theory. In addition to the
relevant local observables like charge and electric field, we compute the
(bipartite) entanglement entropy and show that subtraction of its vacuum value
results in a UV-finite quantity. We find that both string formation and string
breaking leave a clear imprint on the resulting entropy profile. Finally, we
also study the case of fractional probe charges, simulating for the first time
the phenomenon of partial string breaking.
