DEy := (y - 1)^2*(8*y^2*z + y^2 + 8*z^2 - 6*z - 2)*(16*y^2*z - 16*z^2 - 8*z - 1)*(y + 1)^2*y*diff(F(y), y,y,y,y)
  + (y^2 - 1)*(1536*y^6*z^2 + 192*y^6*z + 896*y^4*z^3 - 2160*y^4*z^2 - 1152*y^2*z^4 - 592*y^4*z - 608*y^2*z^3
    - 7*y^4 + 792*y^2*z^2 - 128*z^4 + 278*y^2*z + 32*z^3 + 15*y^2 + 72*z^2 + 22*z + 2)*diff(F(y), y,y,y)
  + 2*(2368*y^6*z^2 + 296*y^6*z + 2816*y^4*z^3 - 3960*y^4*z^2 - 1216*y^2*z^4 - 1020*y^4*z - 2544*y^2*z^3 - 5*y^4
    + 1920*y^2*z^2 + 192*z^4 + 682*y^2*z - 16*z^3 + 15*y^2 - 136*z^2 - 38*z - 2)*y*diff(F(y), y,y)
  + (3712*y^6*z^2 + 464*y^6*z + 7040*y^4*z^3 - 4320*y^4*z^2 - 1024*y^2*z^4 - 1440*y^4*z - 2688*y^2*z^3
    - 2*y^4 + 1056*y^2*z^2 - 384*z^4 + 436*y^2*z + 32*z^3 + 6*y^2 + 272*z^2 + 76*z + 4)*diff(F(y), y)
  + 48*z*(8*y^2*z + y^2 + 24*z^2 - 6*z - 3)*y^3*F(y);


# Typed from paper to verify paper:
L := 2*(1-y^2)*(2*y^2-4*z-1)^2*((1+4*z)^2-16*z*y^2)*Dy^2
	+ 4*y*(2*y^2-4*z-1)*(32*z*y^4-(1+4*z)*((1+28*z)*y^2-4*z*(3+4*z)))*Dy
	+ 16*(2*y^2-14*z-3)*z*y^4+(1+4*z)*(1+32*z+80*z^2)*y^2-(1+4*z)^2*(1+4*z+8*z^2)
	+ (2*(1+8*z)*y^2-4*(1+4*z)*(1-z))*y*SQ:
sq := sqrt(2*(1-4*z)*((1+4*z)^2-16*z*y^2));
L_y_plus  := subs(SQ =  sq, L);
L_y_minus := subs(SQ = -sq, L);

# Verify L_y_plus and L_y_minus:
	with(DEtools): _Envdiffopdomain := [Dy, y]:
	L4 := symmetric_product(L_y_plus, L_y_minus);
	# Check that this is the same as DEy:
	ShouldBeZero := rem(L4, de2diffop(DEy, F(y)), Dy);

chvar := proc(L, y, Dy, sb1,sb2)
	local i, f, x, Dx, s;
	Dx, x := op(_Envdiffopdomain);
	f := add( subs(sb1,sb2,coeff(L, Dy, i)) * mult( Dx/diff(subs(sb1,sb2,y),x) $ i), i = 0 .. degree(L, Dy));
	s := proc(a) factor(simplify(normal(a,expanded),symbolic)) end:
	sort(collect(f/lcoeff(f, Dx), Dx, s), Dx);
end:

from_yz_to_xt :=
	y = (x + (t^2-1)^2) * sqrt(z/x),
	z = 1/(4*(1-2*t^2));

from_xt_to_yz := convert( solve({from_yz_to_xt}, {x,t}), radical);

Divide_by := ( x * (x + (t-1)*(t+1)^3)
                 / (x + (t+1)*(t-1)^3) )^(1/4);

# Compared it with the paper, and it looks identical:
L_plus_x := Dx^2 + ( 1/x - 1/(x + (t+1)*(t-1)^3) + 2*(x+t^4+2*t^2-1)/(x^2 + (2*t^4+4*t^2-2)*x + (t^2-1)^4) )*Dx
  + ( 4*x^2 + (8*t^4-16*t^3-4*t^2+20*t-9)*x+(4*t^4-8*t^3+12*t^2+4*t-5)*(t+1)*(t-1)^3 )
  / ( 16 * x * (x^2+(2*t^4+4*t^2-2)*x+(t^2-1)^4) * (x+(t+1)*(t-1)^3) ) ;
L_minus_x := subs(t = -t, L_plus_x);

# Verify that the change of variables and division sends L_y_plus to L_plus_x:
	_Envdiffopdomain := [Dx,x];
	L := chvar(L_y_plus, y, Dy, from_yz_to_xt):
	L1 := Dx + evala( diff(Divide_by, x)/Divide_by  ):
	ShouldBeZero := normal( symmetric_product(L, L1) - L_plus_x);

L_1_x := Dx - 1/(2*x);
symmetric_product(L_1_x, symmetric_product(L_plus_x, L_minus_x)):
# Check that this matches DEy as claimed in Section 2.2:
chvar(de2diffop(DEy, F(y), [Dy,y]), y, Dy, from_yz_to_xt):
ShouldBeZero := normal(%% - %);