In short, there is no simple formula for the $q$-expansion of the transform of a modular form under an Atkin-Lehner operator $W_Q$, in terms of the $q$-expansion of the original modular form.
The action of $W_Q$ on modular symbols is simply $\{\alpha,\beta\} \mapsto \{W_Q \alpha, W_Q \beta\}$. This is easy to implement, and gives the matrix of $W_Q$ in the basis of homology given by (linear combinations of) Manin symbols. If you are interested in the action on a given newform, you just restrict this matrix to the corresponding Hecke module. Otherwise, you compute a linear combination of such transforms.
This works well for $\Gamma_0(N)$ (trivial character), but not in general, because $W_Q$ doesn't preserve the character.
There is an implementation in Magma, as you have mentioned, with the restriction to cusp forms with trivial character and integral weight $\geq 2$. In Pari/GP there is a completely general implementation, for modular forms of arbitrary weight (including weight 1 and half-integral weight), level and character. In the difficult cases, this relies on a theorem of Borisov and Gunnells, which asserts that a cusp form of weight $>2$ on $\Gamma_1(N)$ can be expressed as a linear combination of pairwise products of Eisenstein series, plus an Eisenstein series. The idea being that the action of $W_Q$ on Eisenstein series is very easy to write down. It should be said, however, that finding such a linear combination of pairwise products, if done exactly, turns out to be very costly. For this reason, the Pari/GP algorithm is numerical; but it works well in practice. Some more details are given in the course Modular forms in Pari/GP by Belabas and Cohen.
There are also theoretical results. In your case the character is trivial, so that newforms are eigenvectors for Atkin-Lehner operators. If the newform $f$ satisfies $a_q(f) \neq 0$ for a prime $q$ dividing the level, then there is a formula for the local eigenvalue at $q$, see for example Atkin, Li, Twists of newforms and pseudo-eigenvalues of $W$-operators, Inv. Math. (1978) 48: 221-244. General theoretical formulas are difficult to get. They seem to be more conveniently stated and proved using the adelic language and the Whittaker model of the newform. The resulting formulas involve $\mathrm{GL}_2$ Gauss sums. For trivial character, there is a completely explicit formula in Nelson, Pitale, Saha, Bounds for Rankin-Selberg integrals and quantum unique ergodicity for powerful levels, J. Amer. Math. Soc. (2014), 27: 147-191, see p. 177. I should say here that $\begin{pmatrix} 0 & 1 \\ -2 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} 2 & 0 \\ 0 & 1 \end{pmatrix}$, so that the problem is reduced to the action of $\sigma = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$. I guess that in your case, the formula will be simpler.
