with(DEtools): _Envdiffopdomain:=[D,x]: read("2.03_transformations.mpl"); read("algo.mpl"); L1 := D^2+x; #L2 := gauge(L1,x,0); L2 := D^2-2*D/x + (2+x^3)/x^2; equiv(L1,L2); G:=r1(x)*D+r0(x); R:=numer(op(2,rightdivision(mult(L2,G),L1))): R:=subs({diff(r0(x),x)=r3(x), diff(r1(x),x)=r4(x)},R): collect(R,D,normal); cc := [coeffs(collect(R,D),D), diff(r1(x),x)-r4(x),diff(r0(x),x)-r3(x)]; cc := [x r1(x) - 2 x r3(x) + 2 r0(x) + x |-- r3(x)| - 2 x r4(x), cc:=ListTools[Reverse](cc): Y:=[seq(op(map(y->`if`(has(y,diff),op(1,y),NULL),indets(op(i,cc)))),i=1..4)]; for k to 4 do ccc := coeff(op(k,cc),diff(op(k,Y),x),1); row[k] := [seq(normal(-coeff(subs(diff(op(k,Y),x)=0,op(k,cc)), op(i,Y),1)/ccc),i=1..4)]; od: A := matrix(4,4,[[op(row[1])],[op(row[2])],[op(row[3])],[op(row[4])]]); V[0] := [1,0,0,0]: for i to 4 do V[i] := map(normal, map(diff,V[i-1],x) - convert(evalm( A &* V[i-1]),list)); od: expand(add(a[i]*V[i],i=0..4)): solve({op(%) , a[4]=1}, {a[0],a[1],a[2],a[3],a[4]}): L := eval(add(a[i]*D^i, i=0..4), %); if indets(L) intersect {a[0],a[1],a[2],a[3],a[4]} <{} then error "The vector V[0] is not cyclic, choose another V[0] and re-run it" fi: La := adjoint(L); c[3] := op(1, expsols( diffop2de(La,y(x)), y(x))); R := collect(c[3]*op(1,leftdivision(L, adjoint(D - normal(diff(c[3],x)/c[3])))), D,normal); Y := normal(expand(add(coeff(R, D, i) * V[i], i=0..3))); map(normal, evalm( map(diff,Y,x) - A &* Y ));