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)): Alg := Ore_algebra:-diff_algebra([dx, x], [ds, s]): 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; phi1 := -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); normal( Ore_algebra:-applyopr(P1, F, Alg) - diff(phi1*F, t) ); P2 := -2*(-19*s^2-9*x+13*s^3+7*s-16*x*s^2+24*x*s) * dx + disc * dx^2; phi2 := -t*(3*s-2)*(s-1)*(2*s^2-4*s*t-s+3*t)*(s-t)^2/(t-x)/q1; normal( Ore_algebra:-applyopr(P2, F, Alg) - diff(phi2*F, t) ); MOrd := Groebner:-MonomialOrder(Alg, plex(dx, ds)): P2_prime := collect(Groebner:-Reduce(P2, {P1}, MOrd), ds, factor); collect(P2_prime, ds, p -> bideg(degree(p, s), degree(p, x)));