f := (1-t)*(1-s)*(1-u)/(1-2*(s+t+u)+3*(s*t+t*u+u*s)-4*s*t*u):
F := normal(eval(f, {t=t/s, u=x/t})/s/t):
q := denom(F):
q1 := q/s/t:
disc := factor(discrim(q1, t)):

# Same as input to RectangularSystem.mpl and SecondIteration.mpl.
P1 :=
  2*(s-1)*(3*s^2-6*s+2) +
  (6*x*s^3-2*s^3-10*x*s^2+s^2-4*x^2*s+10*x*s+3*x^2-4*x) * dx +
  2*s*(3*s-2)*(s-1)^2 * ds:
P2 :=
  -2*(-19*s^2-9*x+13*s^3+7*s-16*x*s^2+24*x*s) * dx +
  disc * dx^2:
# Same as computed by SecondIteration.mpl.
P :=
  4*(576*x^3-801*x^2-108*x+74) * dx +
  (4608*x^4-6372*x^3+813*x^2+514*x-4) * dx^2 +
  (1152*x^5-1746*x^4+475*x^3+121*x^2-2*x) * dx^3:
gamma_ :=
  -243*x^3*(36*x^2-5*x+4) +
  36*(1296*x^4-159*x^3+542*x^2+45*x-2)*x * s +
  (-19730*x^2+56+1008*x+15516*x^4-93084*x^3-93312*x^5) * s^2 +
  (82024*x^2-616-4228*x-42420*x^4+196372*x^3+82944*x^5) * s^3 +
  (-165704*x^2+6924*x+83632*x^4-202768*x^3-27648*x^5+2555) * s^4 +
  (81872*x^3+178716*x^2-2624*x-5112-84096*x^4) * s^5 +
  (17064*x^3-99852*x^2-3952*x+5008+32256*x^4) * s^6 +
  (-18432*x^3+23040*x^2+3048*x-1944) * s^7:

psi[1] := -t*(s-1)*(-6*x*s^2+6*s^2*t-s^2+11*x*s-9*s*t-4*x-x*t+4*t)/(t-x):
psi[2] := -t*(3*s-2)*(s-1)*(2*s^2-4*s*t-s+3*t)*(s-t)^2/(t-x)/q1:

phi[0] := s*(108*x+53)*(3*s-2) / (s-1):
phi[1] := gamma_/2/(-1+s)^2/(24*x*s-16*x*s^2-9*x+s-4*s^2+4*s^3)/(s-x):

Alg := Ore_algebra:-diff_algebra([dt, t], [ds, s], [dx, x]):
MOrd := Groebner:-MonomialOrder(Alg, plex(ds, dx, dt)):
Q := phi[0] + phi[1] * dx:

#`Mgfun/NOT_BUGGED_Reduce` works around a bug in Groebner:-Reduce for
# Maple13, and can be replaced by the latter in Maple14.
`Mgfun/NOT_BUGGED_Reduce`(
  P - Ore_algebra:-skew_product(ds, Q, Alg), [P1, P2], MOrd, 'r', 'q');

A1 := add(factor(normal(coeff(q[1], dx, i)/r))*dx^i, i=0..degree(q[1], dx)):
A2 := add(factor(normal(coeff(q[2], dx, i)/r))*dx^i, i=0..degree(q[2], dx)):

gamma_1 := numer(coeff(A1, dx, 1)):
degree(gamma_1, x), degree(gamma_1, s);

gamma_2 := numer(coeff(A2, dx, 0)):
degree(gamma_2, x), degree(gamma_2, s);

S := factor(normal(Ore_algebra:-applyopr(Q, F, Alg))):
T := factor(normal(Ore_algebra:-applyopr(
  Ore_algebra:-skew_product(A1, psi[1], Alg) +
  Ore_algebra:-skew_product(A2, psi[2], Alg),
  F, Alg))):

denom(S) / (2*s*t*q1^2*disc);
denom(T) / (2*s^2*q1^3*disc^2);

U := numer(S)/(s-t):
seq(degree(U, v), v in [x, s, t]);

V := numer(T)/(s-t):
seq(degree(V, v), v in [x, s, t]);

normal(Ore_algebra:-applyopr(P, F, Alg) - diff(S, s) - diff(T, t));