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.
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:

p0 := -eval(P1, {ds=0, dx=0})/coeff(P1, ds, 1):
p1 := -coeff(P1, dx, 1)/coeff(P1, ds, 1):
q := -coeff(P2, dx, 1)/coeff(P2, dx, 2):

Alg := Ore_algebra:-diff_algebra([dx, x], [ds, s]):
with(Groebner):
MOrd := MonomialOrder(Alg, plex(ds, dx)):
GB := Basis([P1, P2], MOrd):

Q := phi[0](s) + phi[1](s)*dx:
ds_Q := collect(diff(Q, s) + Q*ds, {dx, ds}, distributed);

d := 3:
add(eta[i] * NormalForm(dx^i, GB, MOrd), i=0..d) -
( phi[1](s) * NormalForm(dx*ds, GB, MOrd) +
  phi[0](s) * NormalForm(ds, GB, MOrd) +
  eval(ds_Q, ds = 0) ):
zz := collect(%, dx, normal);

M := 2:
N := 7:
eval(zz,
  { phi[0] = unapply( add(w0[i]*s^i, i=0..M)/(s-1), s ),
    phi[1] = unapply( add(w1[i]*s^i, i=0..N)/2/(s-1)^2/disc, s ) }):
collect(%, dx, numer@normal):
collect(%, {dx, s}, distributed):
sys := {coeffs(%, {dx, s})}:
sol := SolveTools:-Linear(sys,
  {seq(eta[i], i=0..d), seq(w0[i], i=0..M), seq(w1[i], i=0..N)}):

# This complicated-looking final normalisation is session independent.
eval([[seq(eta[i], i=0..d)], [seq(w0[i], i=0..M)], [seq(w1[i], i=0..N)]], sol);
%[1][d+1]:
factor( eval(%%, op(indets(%) minus {x}) = denom(%)) ):
sign(%[1][d+1]):
SOL := factor( expand(% * %%) );

# phi[0]:
phi[0]:= factor( add(SOL[2][i+1]*s^i, i=0..M) ) / (s-1);

# phi[1]:
gamma_ := add(SOL[3][i+1]*s^i, i=0..N);
degree(gamma_, x), degree(gamma_, s);
phi[1] := gamma_ / (2*(s-1)^2*disc):

# Check that P - ds Q is 0 modulo the ideal I = <P1, P2>.
P := add(SOL[1][i]*dx^(i-1), i=1..nops(SOL[1]));
Q := phi[0] + phi[1]*dx:
Groebner:-Reduce(P - Ore_algebra:-skew_product(ds, Q, Alg), [P1, P2], MOrd);