The Magma code below finds all values of tau, excluding those that satisfy
tau*(tau^2-1)*(9*tau^2-1)*(tau^2+6*tau+1) = 0,
for which Z(mu(tau)) has a singularity at some other point besides the known A5-singularity.

STEP 1:

// Gives all (non-excluded) values of tau where mu(tau)
// is singular away from P_0, P_1, and R.

R<u,v,t,Z,x,y,tau> := PolynomialRing(Rationals(),7);
c := -tau^2;
alpha := 0;
Q := -16*c^3;
inv := (x-alpha)*t-1;
inv2 := (x-4*tau^3)*u-1;
inv3 := tau*(tau^2-1)*(9*tau^2-1)*(tau^2+6*tau+1)*v-1;
curve := y^2-(x^3+(1+6*c-3*c^2)*x^2+(-16*c^3)*x);
psi0 := (x-alpha)*Z^6+Q*t;
psi1 := (x-alpha)^2*Z^6+Q^2*t^2;
psi2 := (x-alpha)*Z^5+Q*t*Z;
psi3 := y*Z^5-Q*y*t^2*Z;
psi4 := Z^4+Z^2;
psi5 := (x-alpha)*Z^4+Q*t*Z^2;
psi6 := Z^3;
quant1 := 8*(1-c)*c;
quant2 := 1;
quant3 := 4*c^2*(1 - 4*c + c^2);
f1 := psi0+psi2;
f2 := psi0-psi5;
f3 := psi4+2*psi6;
f4 := quant1*psi0+quant2*psi1+quant3*(psi4-2*psi6);
b1 := 4*(tau-1)^2*tau^2*(1+14*tau+34*tau^2+14*tau^3+tau^4);
b2 := 8*tau^3*(3+28*tau+34*tau^2+28*tau^3+3*tau^4);
b3 := -8*tau^6*(1-52*tau-90*tau^2-52*tau^3+tau^4);
b4 := (1+6*tau+tau^2)^2;
comb := b1*f1+b2*f2+b3*f3+b4*f4;
I := ideal<R|comb,curve,inv,inv2,inv3>;
m := Submatrix(JacobianMatrix([comb,curve,inv]),1,3,3,4);
sing := ideal<R|I,Minors(m,3)>;
time EliminationIdeal(sing,6);

Output is:

Ideal of Polynomial ring of rank 7 over Rational Field
Order: Lexicographical
Variables: u, v, t, Z, x, y, tau
Basis:
[
    tau^11 + 31*tau^10 + 717/2*tau^9 + 2038*tau^8 + 48479/8*tau^7 + 
        73647/8*tau^6 + 53611/8*tau^5 + 20851/8*tau^4 + 4565/8*tau^3 + 
        565/8*tau^2 + 37/8*tau + 1/8
]
Time: 22.460


If one reduces this polynomial modulo p and finds a root, one generally finds that the corresponding curve has 2 additional A_1 singularities, at points above the "other" 2-torsion points of E.  Thus the curve has two A_5 and two A_1 singularities.


STEP 2:

// Gives all (non-excluded) values of tau where Z(mu(tau))
// has a singularity above R that satisfies Z != 1.

// p := 113;
// R<u,v,t,Z,x,y,tau> := PolynomialRing(GF(p),7);
R<u,v,t,Z,x,y,tau> := PolynomialRing(Rationals(),7);
c := -tau^2;
alpha := 0;
Q := -16*c^3;
inv := (x-alpha)*t-1;
inv2 := (Z-1)*u-1;
inv3 := tau*(tau^2-1)*(9*tau^2-1)*(tau^2+6*tau+1)*v-1;
curve := y^2-(x^3+(1+6*c-3*c^2)*x^2+(-16*c^3)*x);
psi0 := (x-alpha)*Z^6+Q*t;
psi1 := (x-alpha)^2*Z^6+Q^2*t^2;
psi2 := (x-alpha)*Z^5+Q*t*Z;
psi3 := y*Z^5-Q*y*t^2*Z;
psi4 := Z^4+Z^2;
psi5 := (x-alpha)*Z^4+Q*t*Z^2;
psi6 := Z^3;
quant1 := 8*(1-c)*c;
quant2 := 1;
quant3 := 4*c^2*(1 - 4*c + c^2);
f1 := psi0+psi2;
f2 := psi0-psi5;
f3 := psi4+2*psi6;
f4 := quant1*psi0+quant2*psi1+quant3*(psi4-2*psi6);
b1 := 4*(tau-1)^2*tau^2*(1+14*tau+34*tau^2+14*tau^3+tau^4);
b2 := 8*tau^3*(3+28*tau+34*tau^2+28*tau^3+3*tau^4);
b3 := -8*tau^6*(1-52*tau-90*tau^2-52*tau^3+tau^4);
b4 := (1+6*tau+tau^2)^2;
comb := b1*f1+b2*f2+b3*f3+b4*f4;
I := ideal<R|comb,curve,inv,inv2,inv3>;
m := Submatrix(JacobianMatrix([comb,curve,inv]),1,3,3,4);
sing := ideal<R|I,Minors(m,3),x-4*tau^3>;
time EliminationIdeal(sing,6);

Output is:

Ideal of Polynomial ring of rank 7 over Rational Field
Order: Lexicographical
Variables: u, v, t, Z, x, y, tau
Homogeneous
Basis:
[
    1
]
Time: 0.120

This means that there are no values of tau that satisfy this.  (Another way to verify this: Remove "inv2" from the definition of I, then compute RadicalDecomposition(sing).  The only things appearing will satisfy Z=1.)


STEP 3:

// Gives all values of tau where Z(mu(tau)) has a
// singularity above P_0, using one coordinate chart.


R<t,V,tau> := PolynomialRing(Rationals(),3);
c := -tau^2;
A := -16*c^3;
psi0 := A;
psi1 := V^6;
psi2 := t*(V*(V^4+A));
psi3 := V*(V^4-A);
psi4 := V^2;
psi5 := V^4;
psi6 := t*V^3;
quant1 := 8*(1-c)*c;
quant2 := 1;
quant3 := 4*c^2*(1 - 4*c + c^2);
f1 := psi0+psi2;
f2 := psi0-psi5;
f3 := psi4+2*psi6;
f4 := quant1*psi0+quant2*psi1+quant3*(psi4-2*psi6);
b1 := 4*(tau-1)^2*tau^2*(1+14*tau+34*tau^2+14*tau^3+tau^4);
b2 := 8*tau^3*(3+28*tau+34*tau^2+28*tau^3+3*tau^4);
b3 := -8*tau^6*(1-52*tau-90*tau^2-52*tau^3+tau^4);
b4 := (1+6*tau+tau^2)^2;
comb := b1*f1+b2*f2+b3*f3+b4*f4;
m := Submatrix(JacobianMatrix([comb]),1,1,1,2);
sing := ideal<R|comb,Minors(m,1)>;
EliminationIdeal(sing,2);

Output is:
[
    tau^32 + 38*tau^31 + 615*tau^30 + 5476*tau^29 + 28691*tau^28 + 
        85422*tau^27 + 109429*tau^26 - 99664*tau^25 - 506598*tau^24 - 
        458708*tau^23 + 367862*tau^22 + 934872*tau^21 + 367862*tau^20 - 
        458708*tau^19 - 506598*tau^18 - 99664*tau^17 + 109429*tau^16 + 
        85422*tau^15 + 28691*tau^14 + 5476*tau^13 + 615*tau^12 + 
        38*tau^11 + tau^10
]

This factors as (-1 + tau)^4 tau^10 (1 + tau)^6 (1 + 6 tau + tau^2)^6, which are the excluded values of tau.

STEP 3:

// Gives all values of tau where Z(mu(tau)) has a
// singularity above P_0, using the other coordinate chart.

R<t,U,tau> := PolynomialRing(Rationals(),3);
c := -tau^2;
A := -16*c^3;
psi0 := A*U^6;
psi1 := 1;
psi2 := t*(U*(1+A*U^4));
psi3 := U*(1-A*U^4);
psi4 := U^4;
psi5 := U^2;
psi6 := t*U^3;
quant1 := 8*(1-c)*c;
quant2 := 1;
quant3 := 4*c^2*(1 - 4*c + c^2);
f1 := psi0+psi2;
f2 := psi0-psi5;
f3 := psi4+2*psi6;
f4 := quant1*psi0+quant2*psi1+quant3*(psi4-2*psi6);
b1 := 4*(tau-1)^2*tau^2*(1+14*tau+34*tau^2+14*tau^3+tau^4);
b2 := 8*tau^3*(3+28*tau+34*tau^2+28*tau^3+3*tau^4);
b3 := -8*tau^6*(1-52*tau-90*tau^2-52*tau^3+tau^4);
b4 := (1+6*tau+tau^2)^2;
comb := b1*f1+b2*f2+b3*f3+b4*f4;
m := Submatrix(JacobianMatrix([comb]),1,1,1,2);
sing := ideal<R|comb,Minors(m,1)>;
EliminationIdeal(sing,2);

Output is:

Ideal of Polynomial ring of rank 3 over Rational Field
Order: Lexicographical
Variables: t, U, tau
Basis:
[
    tau^20 + 36*tau^19 + 542*tau^18 + 4356*tau^17 + 19437*tau^16 + 
        42192*tau^15 + 5608*tau^14 - 153072*tau^13 - 206062*tau^12 + 
        106488*tau^11 + 360948*tau^10 + 106488*tau^9 - 206062*tau^8 - 
        153072*tau^7 + 5608*tau^6 + 42192*tau^5 + 19437*tau^4 + 
        4356*tau^3 + 542*tau^2 + 36*tau + 1
]

This factors as (-1 + tau)^4 (1 + tau)^4 (1 + 6 tau + tau^2)^6.