Fourier expansion at inequivalent cusps

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.

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).

Since this has not yet been mentioned explicitly, I would like to add that the adelic language is a powerful tool for computing Fourier expansions at arbitrary cusps.

Specifically, if $f$ is a cuspidal holomorphic newform then the Fourier coefficients of $f$ at a given cusp are given by values of the Whittaker newform of $f$ at certain adelic matrices. The Whittaker newform associated to $f$ is a global object which factors as a product of local newforms. Concretely, the local newforms are functions on $\mathrm{GL}_2(\mathbb{Q}_p)$ depending only on the local automorphic representation associated to $f$. Since Loeffler and Weinstein have given an algorithm (implemented e.g. in Sage) to compute these local representations, it should be possible to compute the local newforms explicitly. For more details, you can look at the recent article of Corbett and Saha, On the order of vanishing of newforms at the cusps, especially Section 3.

EDIT. Hao Chen has also developed some algorithms to solve this problem in his PhD thesis, Chapter 4.