Let $\Gamma\subset SL(2,\mathbb{R})$ be a Fuchsian group of the first kind. Let $c_1, c_2$ be inequivalent cusps of $\Gamma.$ Consider $f\in M_k(\Gamma)$ a weight $k$ holomorphic automorphic form, and suppose the Fourier expansion of $f$ at the cusp $c_1$ is known.
Given the above expansion, is there an algorithm to compute (even numerically) the Fourier expansion of $f$ at the cusp $c_2$? Thanks!
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).
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.
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.
A chunk of the Princeton PhD thesis of Dan Collins is devoted to this problem. See Collins' preprint here: https://arxiv.org/abs/1802.09740, and the accompanying Sage code.
