  We present an algorithm to evaluate Matsubara sums for Feynman diagrams
comprised of bare Green's functions with single-band dispersions with local U
Hubbard interaction vertices. The algorithm provides an exact construction of
the analytic result for the frequency integrals of a diagram that can then be
evaluated for all parameters $U$, temperature $T$, chemical potential $\mu$,
external frequencies and internal/external momenta. This method allows for
symbolic analytic continuation of results to the real frequency axis, avoiding
any ill-posed numerical procedure. When combined with diagrammatic Monte-Carlo,
this method can be used to simultaneously evaluate diagrams throughout the
entire $T-U-\mu$ phase space of Hubbard-like models at minimal computational
expense.
