# Same as in the previous file:
DEz := (4*z - 1)*z^2*(16*y^2*z - 16*z^2 - 8*z - 1)*(8*z - 1)*diff(F(z), z,z,z,z)
  + z*(4608*y^2*z^3 - 1472*y^2*z^2 - 5888*z^4 + 120*y^2*z - 352*z^3 + 368*z^2 + 2*z - 5)*diff(F(z), z,z,z)
  + (9728*y^2*z^3 - 2624*y^2*z^2 - 17280*z^4 + 196*y^2*z + 144*z^3 + 668*z^2 - 38*z - 4)*diff(F(z), z,z)
  + (4224*y^2*z^2 - 1008*y^2*z - 12480*z^3 + 54*y^2 + 1032*z^2 + 276*z - 18)*diff(F(z), z)
  + (96*y^2*z - 36*y^2 - 960*z^2 + 168*z + 12)*F(z);

# Typed from paper (to make sure the paper doesn't have a typo in the formula):
L := z*(1-4*z)*(1+4*z)^2*((1+4*z)^2-16*y^2*z)*Dz^2
   + (8*z*(48*z^2+8*z-3)*y^2+(1-32*z^2)*(1+4*z)^2)*(1+4*z)*Dz
   + (64*z^3+80*z^2+4*z-1)*y^2 - 3*z*(1+4*z)^3 + (1-8*z)*y * SQ:
sq := sqrt(2*(4*z-1)*(16*z*y^2-(4*z+1)^2));
L_z_plus  := subs(SQ =  sq, L);
L_z_minus := subs(SQ = -sq, L);

with(DEtools): _Envdiffopdomain := [Dz,z]:
L4 := symmetric_product(L_z_plus, L_z_minus);
# Check that this is the same as DEz:
ShouldBeZero := rem(L4, de2diffop(DEz, F(z)), Dz);