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.

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.