Page 523
# first define the old variables
p123:=(a0*b0*c0*d0+3*a0*b1*c1*d1+
%ad 3*a1*(b1*c0*d0+b0*c1*d1+2*b1*c1*d1));
p12:=(12/4)*(a0*(b0*c0*d1+b1*c1*d0+2*b1*c1*d1)+ a1*(b0*c1*d0+
 b1*c0*d1+2*b1*c1*d1)+2*a1*(b1*c0*d1+b1*c1*d0+b0*c1*d1+b1*c1*d1));
p13:=(12/4)*(a0*(b0*c1*d0+b1*c0*d1+2*b1*c1*d1)+ a1*(b1*c1*d0+
 b0*c0*d1+2*b1*c1*d1)+2*a1*(b1*c1*d0+b1*c0*d1+b0*c1*d1+b1*c1*d1));
p23:=(12/4)*(a1*(b0*c0*d0+3*b1*c1*d1)+ a0*(b1*c0*d0+b0*c1*d1+
 2*b1*c1*d1) +2*a1*(b1*c0*d0+2*b1*c1*d1+b0*c1*d1));
pdis:=(24/4)*(a0*(b0*c1*d1+b1*c0*d1+b1*c1*d0+b1*c1*d1)+
a1*(2*b1*c1*d1+b0*c0*d1+b1*c1*d0)+ a1*(2*b1*c1*d1+b1*c0*d1+b0*c1*d0)
 +a1*(b1*c1*d1+b1*c0*d1+b1*c1*d0+b0*c1*d1));

# now compute the new variables
P123:=algsubs(a0*b0=e0-3*a1*b1,p123): P123:=expand(P123):
 P123:=algsubs(a0*b1=e1-a1*b0-2*a1*b1,P123);
P12:=expand(p12): P12:=algsubs(a0*b0=e0-3*a1*b1,P12):
 P12:=expand(P12): P12:=algsubs(a0*b1=e1-a1*b0-2*a1*b1,P12);
P13:=expand(p13): P13:=algsubs(a0*b0=e0-3*a1*b1,P13):
 P13:=expand(P13): P13:=algsubs(a0*b1=e1-a1*b0-2*a1*b1,P13);
P23:=expand(p23): P23:=algsubs(a0*b0=e0-3*a1*b1,P23):
 P23:=expand(P23): P23:=algsubs(a0*b1=e1-a1*b0-2*a1*b1,P23);
Pdis:=expand(pdis): Pdis:=algsubs(a0*b0=e0-3*a1*b1,Pdis):
 Pdis:=expand(Pdis): Pdis:=algsubs(a0*b1=e1-a1*b0-2*a1*b1,Pdis);

              Page 526
f1:=(1-3*t1)*(1-3*t2)*(1-3*t3)+3*t1*t2*t3;
f2:=3*(1-3*t1)*(1-3*t2)*t3+3*t1*t2*(1-3*t3)+6*t1*t2*t3;
f3:=3*(1-3*t1)*t2*(1-3*t3)+3*t1*(1-3*t2)*t3+6*t1*t2*t3;
f4:=3*(1-3*t1)*t2*t3+3*t1*(1-3*t2)*(1-3*t3)+6*t1*t2*t3;
f5:=6*(1-3*t1)*t2*t3+6*t1*(1-3*t2)*t3+6*t1*t2*(1-3*t3)+6*t1*t2*t3;
cden:=f1*f2*f3*f4*f5;
# fix t3 to reduce computation time, then p3 won't be needed
t3:=1/10;
u1:=31: u2:=5: u3:=7: u4:=11: u5:=13:
p1:=u1*diff(f1,t1)*cden/f1 + u2*diff(f2,t1)*cden/f2 +
  u3*diff(f3,t1)*cden/f3 + u4*diff(f4,t1)*cden/f4 +
  u5*diff(f5,t1)*cden/f5;
p2:=u1*diff(f1,t2)*cden/f1 + u2*diff(f2,t2)*cden/f2 +
  u3*diff(f3,t2)*cden/f3 + u4*diff(f4,t2)*cden/f4 +
  u5*diff(f5,t2)*cden/f5;
# the next for the Hessian matrix
h11:=cden*diff(p1,t1)-p1*diff(cden,t1);
h12:=cden*diff(p1,t2)-p1*diff(cden,t2);
h21:=cden*diff(p2,t1)-p2*diff(cden,t1);
h22:=cden*diff(p2,t2)-p2*diff(cden,t2);
simplify(h12-h21); # check that the mixed partials are equal
# use fsolve to find a root
S1:=fsolve({ p1,p2},{ t1,t2},{ t1=0..1,t2=0..1});
# see if root of cden, check size of cden against p1 and p2
assign(S1); p1; p2; cden;
# if this is a root of cden, try another
t1:='t1'; t2:='t2'; # reset t1 and t2
S2:=fsolve({ p1,p2},{ t1,t2},{ t1=0..1,t2=0..1},avoid={ S1});
# and avoid={ S1,S2} if another round necessary
# e.t.c. also check eigenvalues of hessian are negative for maximum
assign(S2);
h:=array([[h11,h12],[h21,h22]]);
evalf(Eigenvals(h));