Suppose $f$ is a classical modular form of weight $r$ for a (congruence) group $\Gamma$. Let $\gamma$ be any matrix in $\operatorname{SL}\_2(\mathbb{Z})$. Then the slash operator $|\_\gamma$ is usually defined as:
$$f|\_\gamma(z) = \frac{f(\gamma z)}{(cz+d)^r}$$
where $(c, d)$ stands for the bottom row of $\gamma$. The function $f|\_\gamma$ is again a modular form of the same weight but for the conjugated group $\gamma^{-1} \Gamma \gamma$. I'm interested in some method to compute the first $10^5$ Fourier coefficients of $f|\_\gamma$. The ones for $f$ are known. I'm specially interested in the case where $f$ is the cuspform associated to the elliptic curve $y^2 + yx + y = x^3 + 4x - 6$ (group $\Gamma\_0(14)$).


I do not know if this is a simple task or not, as I barely know anything about computational algorithms for modular forms, but simply approximating the scalar product against an exponential function with good precision is probably a tricky approach.



Computing coefficients for the slash operator of a modular form
If you recall the definition of a modular form, you ask it to be holomorphic at all cusps. Computing the Fourier expansion at the infinity cusp of the form $f|\_\gamma$ ammounts to compute the q-expansion at the cusp $\gamma \cdot \infty$ of the original form. There is a nice trick (due to Asai, <https://projecteuclid.org/euclid.jmsj/1240434239>) which works for any form of level $\Gamma\_0(N)$ if $N$ is square free which is due to the fact that the Atkin-Lehner involutions permute the cusps transitively in this case. You can look at the details in the book "How One Can Compute in Polynomial Time the Value of Ramanujan's Tau at a Prime", by Bas Edixhoven and Jean-Marc Couveignes, page 144 (they present the solution as an actual algorithm). I do not know if this has been implemented in magma or sage.


If $N$ is not square free, you need some numerical method as you mention in general.


