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.