This can be done numerically.  The $n$-th Fourier coefficient around the cusp is given by an integral of the form $$\int_0^{1} f|_{\sigma_2} (z) e(-nz) dx,$$ where $\sigma_2$ is a scaling matrix for the cusp $c_2$.  See p.43 of Iwaniec's Topics in Classical Automorphic Forms for definitions.  In the above, $z=x+iy$ and the integral is independent of $y>0$ which generally one picks to make the calculation as efficient as possible.  Also, I am assuming that the multiplier system is such that $f|_{\sigma_2}$ is periodic with period $1$.

In principle, the Fourier expansion of $f$ allows one to numerically approximate $f|_{\sigma_2}$ at any $z \in \mathbb{H}$, so the above integral can be numerically calculated to any degree of precision.

added later:  The same method works for non-congruence subgroups, as well as more general automorphic forms (e.g. Maass cusp forms).