with(DEtools):_Envdiffopdomain:=[Dx,x];

# Integrate this function
f:=ln(x)^2*x;

f:=erf(x);

f:=exp(x)*sin(x)^2;


# Find operator M for which M(f)=0. Search up to order 6.
for n to 6 do
	M:=add(a[i]*Dx^i,i=0..n);
	Mf:=numer(normal(add(a[i]*diff(f,[x$i]),i=0..n)));
	v:=indets(Mf) minus indets([M,x]);
	while v<>{} do
		Mf:=seq(coeffs(collect(i,v[1]),v[1]),i=[Mf]);
		v:=v minus {v[1]}
	od;
	S:=solve({a[n]-1,Mf}, indets(M) minus {Dx});
	if S<>NULL then
		M:=subs(S,M);
		break
	fi;
od;

if n<=6 then
  # Found operator.
  v:=integrate_sols(M);
  if type(v,list) then
	r:=v[2];
	result:=add(coeff(r,Dx,i)*diff(f,[x$i]),i=0..degree(r,Dx));
	print(result);
  fi;
fi;