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.
ADDED: complement to David Loeffler's answer: the above method needs a number of coefficients proportional to $N$, BUT is valid only for $k\ge2$ integral. Dan Collins' method on the other hand is also valid for $k=1$ and $k$ half-integral, but requires a number of coefficients proportional to $N^2$. It is also available in Pari/GP 2.11.