We developed an analytical theory for the many-body potential of mean force (POMF) between N spheres immersed in a continuum chain fluid. The theory is almost exact for a T polymer solution in the protein limit (small particles, long polymers), where N-body effects are important. Polydispersity in polymer length according to a SchulzFlory distribution emerges naturally from our analysis, as does the transition to the monodisperse limit. The analytical expression for the POMF allows for computer simulations employing the complete N-body potential (i.e., without n-body truncation; n < N). These are compared with simulations of an explicit particle/polymer mixture. We show that the theory produces fluid structure in excellent agreement with the explicit model simulations even when the system is strongly fluctuating, e.g., at or near the spinodal region. We also demonstrate that other commonly used theoretical approaches, such as truncation of the POMF at the pair level or the Asakura Oosawa model, are extremely inaccurate for these systems.