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));