I am really embarrassed to ask this question since I am supposed to be an expert on modular forms, but here goes anyway since I have been stuck on this for several days.

Let $q$ be a prime divisor of $N$ such that $q^2\nmid N$, and let $\chi$ be some Dirichlet character modulo $N$. If $F$ is a normalized eigenform (newform of course) in $S_k(\Gamma_0(N),\chi)$ with Fourier coefficients $a(n)$, it is known from Atkin--Lehner--Li that if $W_q$ denotes an Atkin--Lehner involution then $F$ is an eigenform for the action of $W_q$ with eigenvalue $-a(q)/q^{k/2-1}$, and also that $a(q)^2=\chi(q)q^{k/2-1}$.

I am implementing this in Pari/GP and found a contradiction. Most probably my implementation is wrong, but please follow the following simple example (I know that David Loeffler had a ticket about a similar problem in Sage, but that was solved).

Consider the space $S_3(\Gamma_0(12),\chi_{-4})$, where as usual $\chi_{-4}$ is the odd character mod 4. It has a single eigenfunction defined over $\mathbb Q(\sqrt{-3})$, together with its conjugate. I can choose $W_3=[3,-1;12,-3]$, and $\tau=i=\sqrt{-1}$, so that $W_3(\tau)=(3i-1)/(12i-3)$. I then compute $1000$ Fourier coefficients of $F$, and the value $F(i)$ and $F|_3W_3(i)$. Instead of finding $-a(3)/3^{1/2}$ as ratio, I find $+a(3)/3^{1/2}$ (equal to $i$ and $-i$, but in the wrong order or equivalently with the wrong sign).

I did many other experiments, both in odd and even weights, and with nonreal characters, and my implementation tells me that the eigenvalue should be $-a(q)\overline{\chi(q)}/q^{k/2-1}$.

Now I went through Winnie Li's proof which I essentially reproduced in my book with F. Stromberg, and it seems correct. So I am stuck. Can someone at least check the above example, or tell me what I am doing wrong ? Note that I am not mixing the two possible embeddings of my form: both give me the wrong sign.

ADDED March 27: I am very surprised not to have any (positive or negative) answer. Since the error is systematic, let me give another example in even weight with numerical values. Consider $S_2(\Gamma_0(35),\chi_5)$. It is of dimension $2$ generated by an eigenform and its conjugate. One specific embedding of the eigenform is (sorry for the GP notation) $$F(\tau)=q + 2*I*q^2 - I*q^3 - 2*q^4 + (-2 - I)*q^5 + 2*q^6 - I*q^7 + 2*q^9 + O(q^{10}).$$ I take a random point, say $\tau=I/2+1/Pi$. Simply by summing sufficiently many terms I find $F(\tau)=-0.015178...+0.03676...*I$.

I choose as Atkin--Lehner matrix for the prime $q=7$ the matrix $W_7=[7,-3;35,-14]$ which has the required properties. I compute $F|_2W_7(\tau)$ (here I use 2000 terms since the imaginary part is small) and I find that it is equal to $0.03676+0.015178...*I$, so the eigenvalue is $-I$. However, Atkin--Lehner-Li tells me it should be equal to $-a(7)$, which as one sees from above is equal to $+I$, contradiction. Can someone please check this in Sage, Magma, or for that matter Pari/GP ?