Pari/GP 2.11 does this for $\Gamma_1(N)$ (more precisely for the spaces $M_k(\Gamma_0(N),\chi)$). The algorithm is based on a variant of a theorem of Borisov--Gunnells on the generation by products of two Eisenstein series.