interface(quiet=true): with(linalg):

f := u0 + u1*x + u2*x^2 + u3*x^3 + u4*x^4;

dis := 256*u0^3*u4^3-192*u0^2*u1*u3*u4^2-128*u0^2*u2^2*u4^2+144*u0^2*u2*u3^2*u4-27*u0^2*u3^4+144*u0*u1^2*u2*u4^2-6*u0*u1^2*u3^2*u4-80*u0*u1*u2^2*u3*u4+18*u0*u1*u2*u3^3+16*u0*u2^4*u4-4*u0*u2^3*u3^2-27*u1^4*u4^2+18*u1^3*u2*u3*u4-4*u1^3*u3^3-4*u1^2*u2^3*u4+u1^2*u2^2*u3^2;

A := array([
[a11,a12,a13,a14],
[a12,a22,a23,a24],
[a13,a23,a33,a34],
[a14,a24,a34,a44]]):

B := array([
[b11,b12,b13,b14],
[b12,b22,b23,b24],
[b13,b23,b33,b34],
[b14,b24,b34,b44]]):

h := sort(collect(det(matadd(scalarmul(A,x),B)),[x],distributed)):
U0 := coeff(h,x,0):  U1 := coeff(h,x,1): 
U2 := coeff(h,x,2):  U3 := coeff(h,x,3):  U4 := coeff(h,x,4): 
ss := {u0=U0,u1=U1,u2=U2,u3=U3,u4=U4}:
for s in ss do lprint(s); print(` `); od:
print(nops(U0),nops(U1),nops(U2),nops(U3),nops(U4)); print(` `);

hurwitz := 0:

for i from 1 to nops(dis) do
ergo := sort(expand(subs(ss,op(i,dis)))):
lprint(i,op(i,dis), time(), nops(ergo));
hurwitz := hurwitz + ergo:
od:

degree(hurwitz);
nops(hurwitz);