Here is Magma code that gets you the answer in a few seconds. I made a special case for the bad primes, and did them by hand.
_<x> := PolynomialRing(Rationals());
f5 := 344 + 3106*x - 1795*x^2 - 780*x^3 - x^4 + x^5;
g24 := 14488688572801 - 2922378139308818*x^2 + 134981448876235615*x^4 -
1381768039105642956*x^6 + 4291028045077743465*x^8 -
2050038614413776542*x^10 + 287094814384960835*x^12 -
9040633522810414*x^14 + 63787035668165*x^16 - 158664037068*x^18 +
152929135*x^20 - 50726*x^22 + x^24;
K := NumberField(f5);
_,D := IsSquare(Integers()!Discriminant(f5));
prec := 30;
CHAR_TABLE := CharacterTable(GaloisGroup(g24));
chi := CHAR_TABLE[2];
BAD_FACTORS :=
[ <2,Polynomial([1,-1,1])>,
<3,Polynomial([1,-ComplexField(prec)!chi[9],1])>,
<5,Polynomial([1,0,1])>,
<7,Polynomial([1,0,1])>,
<71,Polynomial([1,0,1])>,
<137,Polynomial([1,1,1])>,
<163,Polynomial([1,-ComplexField(prec)!chi[5],1])>,
<1951,Polynomial([1])>,
<16061,Polynomial([1,-2,1])>,
<889289,Polynomial([1,-ComplexField(prec)!chi[8],1])> ];
BAD := [bf[1] : bf in BAD_FACTORS];
FACTORS := [bf[2] : bf in BAD_FACTORS];
function LOCAL(p,d : Precision:=prec)
if p in BAD then return FACTORS[Position(BAD,p)]; end if;
R := Roots(ChangeRing(f5,GF(p)));
if #R eq 1 then return Polynomial([1,0,1]); end if;
if #R eq 2 then
ord := Lcm([Degree(f[1]) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 3 select 1 else -1,1]); end if;
if #R eq 5 then
ord := Lcm([Degree(f[1]) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 1 select -2 else 2,1]); end if;
r := Roots(ChangeRing(f5,GF(p^5)));
x := r[1][1];
prod := GF(p)!&*[x^(p^i)-x^(p^j) : j in [(i+1)..4], i in [0..4]];
wh := prod eq GF(p)!D;
ord := Lcm([Degree(f[1]) : f in Factorization(Polynomial(GF(p),g24))]);
if ord eq 10 then class := wh select 8 else 9; // compatible with FACTORS
else class := wh select 6 else 5; end if;
return Polynomial([1,-ComplexField(prec)!chi[class],1]);
end function;
L := LSeries(1, [0,0], 1951^2, LOCAL : Precision:=prec);
// s->1-s, Gamma(s/2)^2
psi := DirichletGroup(1951, CyclotomicField(10)).1;
p1951 := Polynomial([1,-ComplexField(prec)!CyclotomicField(5).1]);
TP := TensorProduct(L, LSeries(psi : Precision:=prec), [<1951, 1, p1951>]);
CheckFunctionalEquation(TP);
Here is the special values:
ev := Evaluate(TP,0); // 2-1.453085056...
rel := PowerRelation(ev,4 : Al:="LLL");
NF := NumberField(rel);
Q5<zeta5> := CyclotomicField(5);
assert IsIsomorphic(NF,Q5);
Q5!NF.1;
So $L(\rho,0)=-4\zeta_5(1+\zeta_5)$ for Marty. I get $L(\rho_0,-1)=32(48723\sqrt{5} - 778741)$ as an algebraic. I get $L(\rho,-2)=8800\zeta_5^3 - 14444\zeta_5^2 + 35604\zeta_5 + 17412$ with more precision. I determined the TensorProduct factor at 1951 via trial and error, making the obvious guesses until one worked (the failure is at 100-110 digits). With this, I take it to 240 digits and I can even get $$L(\rho,-4)=-18475535360\zeta_5^3 - 11142861380\zeta_5^2 - 12091894020\zeta_5 - 7107607296$$ and $$L(\rho,-6)=25255057273186244\zeta_5^3 - 1015274469604000\zeta_5^2 - 15695788409197884\zeta_5 +
9459547822189412$$ The precision can go higher if you want more.
Finally, the Maass form:
function MaassEval(L,z)
x:=Real(z); y:=Imaginary(z);
printf "Using %o coefficients\n", Ceiling(11/y);
C := LGetCoefficients(L,Ceiling(11/y));
pi := Pi(RealField());
a := Sqrt(y)*&+[C[n]*KBessel(0,2*pi*n*y)*Sin(2*pi*n*x) : n in [1..#C]];
return a;
end function;
zz:=0.0001+0.0001*ComplexField().1;
MaassEval(TP,zz);
// Using 110000 coefficients
// -1.71477211817772949974178783985E-8 + 9.01673609747756708674470686948E-9*i
MaassEval(TP,zz/(1951*zz+1));
// Using 161297 coefficients
// -1.71477211817772949974179078240E-8 + 9.01673609747756708674496293450E-9*i
Step I: Put the degree 24 polynomial into Magma, make it a number field, and call LSeries on it. This divides the $L$-function into a product of 7 distinct ones (Dokchitsers code, under an attribute called "prod" on the L-series object), given by Artin representations. So my plan was to compute zeros for each of these $L$-functions, with them being a subset of those of $\zeta_M$ naturally, and recombine.
As representations for $SL(2,3)$ this is as
$$1\oplus \omega\oplus\bar\omega\oplus 2\tau_2\oplus 2\tau_2\omega\oplus 2\tau_2\bar\omega\oplus 3\kappa_3.$$
Or as an $L$-function product
$$\zeta_M=\zeta\cdot L(\omega)\cdot L(\bar\omega)\cdot L(\tau_2)^2\cdot
L(\tau_2\omega)^2\cdot L(\tau_2\bar\omega)^2\cdot L(\kappa_3)^3.$$
EDIT: Oh I see now, you wanted it for $M$ w/o assuming analyticity of Artin constituent parts, but Magma automatically dissembles it, and assumes. But now at this point, I see you were trying to avoid this decomposition perhaps, but then I really don't see how you could work with $\zeta_M$ directly, as the conductor is too big, being $163^{16}$. From the standpoint of scientific evidence, it is probably enough that CheckFunctionalEquation in Magma for each of the above constituents gives an answer indistinguishable from zero, though the feel of zeros is also nice.
Another way to test individual analyticity is to use subfields. There is a unique (up to isomorphism) quartic subfield $K_4\subset M$, with
$$\zeta_K=\zeta\cdot L(\kappa_3).$$
The conductor of $\zeta_K$ is small enough to compute with it directly (no decomposition), so the analyticity of $L(\kappa_3)=\zeta_K/\zeta$ can be checked (numerically) by seeing if Riemann zeros are a subset of those of $\zeta_K$. EDIT: Of course, $\kappa_3$ is induced from the a nontrivial linear character of the $Q_8$ subgroup of $SL(2,3)$, so analyticity already follows by induction.
There is also a unique (up to isomorphism) sextic subfield $L_6\subset M$ with
$$\zeta_L=\zeta\cdot L(\omega)\cdot L(\bar\omega)\cdot L(\kappa_3),$$
so the linear parts can be isolated here by $\zeta_L/\zeta_K$, though maybe superfluous ("easy" theorem)? In that matter, already the cubic subfield $C_3\subset M$ has $\zeta_C=\zeta\cdot L(\omega)\cdot L(\bar\omega)$ if desired.
There is a unique (up to isomorphism) octic subfield $N_8\subset M$ with
$$\zeta_N=\zeta\cdot L(\tau_2\omega)\cdot L(\tau_2\bar\omega)\cdot L(\kappa_3).$$
The integer ring has discriminant $163^4$ so $\zeta_N$ is likely still computable directly. By quotient $\zeta_N/\zeta_K$, this can tell about the product of conjugate 2-dimensional $L$-functions. EDIT: see below for another idea, using twists of $\zeta_N$.
Finally the unique duodecic subfield $R_{12}\subset M$ has
$$\zeta_R=\zeta\cdot L(\omega)\cdot L(\bar\omega)\cdot L(\kappa_3)^3.$$
So this gives nothing new, and the discriminant is too large anyway.
Note that none of the parts has $L(\tau_2)$ directly, only $\zeta_M$ itself that is too hard to compute. EDIT: However, by Rankin-Selberg I think(?) it follows that the analyticity of $L(\tau_2)$ is equivalent to that after twisting by $\omega$ to get $L(\tau_2\omega)$.
Answer?: For that matter, recurring to the octic subfield $N_8$, instead of using just the Dedekind $\zeta$-function of $N_8$, twisting it by $\omega$ could be profitable, achieving
$$L(\sigma_8\omega)=L(\omega)\cdot L(\tau_2\bar\omega)\cdot L(\tau_2)\cdot L(\kappa_3),$$
$$L(\sigma_8\bar\omega)=L(\bar\omega)\cdot L(\tau_2)\cdot L(\tau_2\omega)\cdot L(\kappa_3),$$
where $\sigma_8=1\oplus\tau_2\omega\oplus\tau_2\bar\omega\oplus\kappa_3$ is for the Dedekind of $N_8$. Also $\kappa_3\omega=\kappa_3$ nicely. Noting here the twisting gives
$$L(\sigma_8)L(\sigma_8\omega)L(\sigma_8\bar\omega)=\zeta_M,$$
this could provide a way to compute $\zeta_M$, as each part of the left is known as analytic I presume. The conductor of $L(\sigma_8)$ is $163^4$, and that for the twists is $163^6$. Or one can avoid $\zeta_M$ alternatively, as my computation is
$$L(\tau_2)^2={L(\sigma_8\omega)L(\sigma_8\bar\omega)\zeta^2\over\zeta_{N_8}\zeta_{L_6}},$$
where each factor on the right should be known (easily?) to be holomorphic away from the $\zeta$-pole. Note this exploits the linear characters of $SL(2,3)$, and you have none for your next case of $SL(2,5)$.
In all cases, zeros need to be computed, and the right tool is L-calc. But I don't know if it is really feasible to go too far for $L(\sigma_8\omega)$ of conductor $163^6$, without intensive effort.
Part II: Zeros of $L(\tau_2)$ computation (example):
I compute its first few zeros for the 2-dimensional representation of real character. For this representation, with $10^5$ coeffs (taking 6sec in Magma), I exported these to Lcalc (a somewhat hackish tool of M. O. Rubinstein, included in Sage I suppose but I did it direct), which returns after 12 seconds, listing the imaginary parts of the first 10 zeros on the half-line (also proving RH up to that point by A Turing method):
0
0.99014365233
1.38830360231
2.35103235859
3.45296629741
4.32708276131
4.73989005257
5.42392092883
5.99574967707
6.70167394143
The first zero is at $s=1/2$ as the sign is $-1$. The second is about $s\approx 1/2+0.990i$
For reference, here are the L-calc settings I used, hackish as I say:
1
0
100000
0
2
.5
0.5 0.0
.5
0.5 0.0
51.884511447957879460656106859439682023
-1.0 0.0
0
As explained in their help, the first "1" says the coeffs are integers, the second "0" says they are not special, the 3rd "100000" says they are that many, the 4th "0" says not periodic, the 5th "2" says degree 2, the 6th ".5" and "0.5 0.0" say the form of the gamma factor, the "51.884" is $\sqrt{163^2/\pi^2}$ as the analytic conductor, the "-1.0 0.0" is the sign, and the final "0" says no poles. Then the 100000 coeffs are given as integers.
Best Answer
The algorithm to embed function fields, ie. to test if a function field E can be embedded into a function field F has been developed (and implemented in Magma) by a student of Florian Hess: Gerriet Möhlmann as part of his Diploma work. His thesis (in German) can be found at http://www.math.tu-berlin.de/~kant/publications/diplom/moehlmann.pdf. The method is an extension of Florian's automorphism algorithm, in Magma, it is available through the Inclusions command.
To generate function fields of a given genus there are a few possibilities, none of them worked out completely. If the field can be obtained as an Abelian extension (eg. (hyper)elliptic curves have a degree 2 model) then class field theory can be used to generate all such fields. Similar, soluble exxtensions can be constructed this way. For general extensions, one could use Hunter's theorem to get bounds on the valuations of a primitive element and then enumerate all polynomials that might have such roots. Both methods have in common that they produce too many field extensions that correspond to isomorphic curves. The class field theoretic approach has the advantage of being available through Magma.... (I can provide details if anyone is interested)