[these are comments, not an answer, but there were too many for the comments box]
Hey---I wrote that code too! I did it to teach myself "practical Maass forms". I wrote in pari, not sage. I didn't do the example you did. Here's what I did, for what it's worth. First I tried a dihedral example. I used the Hilbert class field of $\mathbf{Q}(\sqrt{145})$; the class group is cyclic of order 4, giving a $D_8$ extension of $\mathbf{Q}$ with a faithful 2-dimensional representation. It's an easy exercise in factoring polynomials mod $p$ to compute traces of Frobenius, and I got the Hecke eigenvalues with little trouble.
But here's the big question: how do you know you got them right? Here's how I did it. I computed the first 200,000 Hecke eigenvalues, created the formal power series defining a function on the upper half plane as per usual, and then I evaluated it to many decimal places at lots of random points $z$ and $\gamma.z$ with $\gamma$ in the level. In all the cases I tried, the answers were the same to within experimental error. I concluded that probably I'd got everything working.
I also did an $S_3$ extension (the class group of $\mathbf{Q}(\sqrt{79})$) and a non-algebraic example coming from a Grossencharacter---this one had conductor 8 but eigenvalue not $1/4$.
I also tried to do an $A_5$ example! As you probably know from the $A_4$ example, an issue that needs resolving here is that the image of Galois in $GL_2(\mathbf{C}$) isn't $A_4$, it's the central extension, so you need to know how primes split in that last extension, which is computationally more expensive. Bjorn Poonen showed me a wonderful trick though, so I could do it. I computed $a_n$ for hundreds and thousands of $n$, built the function on the upper half plane, and checked to see if it was invariant by the level group. It wasn't :-( I concluded that either the Langlands program was wrong or my code was wrong, and I had a good idea which.
Here is the pari script for the 145 example, by the way:
N=200000;
f=x^4 - x^3 - 3*x^2 + x + 1;
ap(p)=if(p==5,-1,if(p==29,-1,if(issquare(Mod(p,5))&&issquare(Mod(p,29)),2*(matsize(factormod(f,p))[1]-3),0)));
chi(n)=kronecker(n,145);
v=vector(N,i,0);
v[1]=1;
for(i=2,N,fac=factor(i);k=matsize(fac)[1];\
if(k>1,v[i]=prod(j=1,k,v[fac[j,1]^fac[j,2]]),\
if(fac[1,2]==1,v[i]=ap(i),\
p=fac[1,1];e=fac[1,2];v[i]=v[p]*v[p^(e-1)]-chi(p)*v[p^(e-2)]))\
);
F(z)=local(x,y,M);x=real(z);y=imag(z);M=ceil(11/y);if(M>N,error("y too small."));sqrt(y)*sum(n=1,M,if(v[n]==0,0,v[n]*besselk(1e-30*I,2*Pi*n*y)*cos(2*Pi*n*x)))
That's it! It's pretty self-explanatory. ap(p) returns the coefficient $a_p$ of the form. chi is the character of the form. If you run this the computer will pause for a few seconds while it computes the first 200,000 coefficients of the Maass form. After that it will give you a function $F$ on the upper half plane, which is defined by a Fourier expansion, and the miracle will be that it will be $\Gamma_1(145)$-invariant. For example, after running the code above, you can try this:
gp > z=-0.007+0.08*I
%5 = -0.007000000000000000000000000000 + 0.08000000000000000000000000000*I
gp > F(z)
%6 = 0.2101332751524672135753981488 + 0.E-30*I
gp > F(z/(145*z+1))
%7 = 0.2101332751524672135753981489 + 0.E-30*I
gp > %6-%7
%8 = -5.67979851 E-29 + 0.E-30*I
What this says is that for $z$ a random element of the upper half plane such that $z$ and $\gamma.z$ both have imaginary part which is not too small (if the im part is too small you need more Fourier coeffts), where here $\gamma=(1,0;145,1)$, $F$ evaluates, to within experimental error, to the same value at $z$ and $\gamma.z$.
Note that Marty's example is of $A_4$ type, so more interesting than this example, but a theorem of Langlands tells us that Marty's example really will be a cusp form.