Algorithm: Detection of symmetric powers of second order operators.

Notation:    Dx = d/dx
Input:       L = Dx^n + An * Dx^(n-1) + Bn * Dx^(n-2) + ...
Output is either:
	"L is not a symmetric power of a 2nd order operator"
or:
	A,B  such that L is (n-1)-th symmetric power of Dx^2 + A*Dx + B

Step 1:  A := 2/(n*An/(n*(n-1))

Step 2:  B := 6/((n-1)*n*(n+1)) * Bn - (n-2)/(n+1) * An' - (3*n-1)*(n-2)/(n^2*(n-1)^2*(n+1)) * An^2

Step 3:  Compute (n-1)'th symmetric power of Dx^2+A*Dx+B with the algorithm of Bronstein/Mulders/Weil.

Step 4:  Test if L equals this symmetric power.
         If so, then output A,B, otherwise L is not a symmetric power of a 2nd order operator.


How the formulas given in Step 1 and Step 2 can be found is shown in the Maple worksheet on:
	http://www.math.fsu.edu/~hoeij/papers/symprod/sympow_comment.html
but they were already known before, see:
	R. Chalkley, Relative invariants for homogeneous linear differential equations,
	J. Differential Equations 80, 107-153, 1989.
or:
	G. Fano, Ueber lineare homogene Differentialgleichungen mit algebraischen
	Relationen zwischen den Fundamentallosungen, Math. Ann., 53, 493-590, 1900.


To see the implementation, type this in Maple 5:
	interface(verboseproc=2):
	readlib(`dsolve/diffeq/higherorder`);

and one sees:

	proc(A, x)
	local i, a, DF, newA, b, ans, m, compare, y;
	option `Copyright (c) 1997 George Labahn. All rights reserved.`;
	    m := nops(A) - 2;
	    newA := [seq(A[i]/A[m + 2], i = 1 .. m + 2)];
	    a := 2*newA[m + 1]/(m*(m + 1));
	    b := normal(1/4*(-3*m^4*a^2 - 2*m^3*a^2 - 4*m^3*diff(a, x) + 3*m^2*a^2
	         + 2*m*a^2 + 4*m*diff(a, x) + 24*newA[m])/(m*(m^2 + 3*m + 2)));
	    compare := DEtools[symmetric_power](DF^2 + a*DF + b, m, [DF, x]);
	    for i from 0 to m do
	        if simplify(coeff(compare, DF, i) - newA[i + 1]) <> 0 then
	            RETURN([])
	        fi
	    od;
	    ans := readlib('`dsolve/diffeq/linearODE`')(b*y + a*_.y[1] + _.y[2], x,
	        y);
	    if ans <> [] and nops(ans) = 2 and not has(ans, 'DESol') then
	        RETURN([seq(ans[1]^i*ans[2]^(m - i), i = 0 .. m)])
	    fi;
	    []
	end

To see the code in more recent Maple versions, type:
	interface(verboseproc=2):
	op(`dsolve/diffeq/higherorder/is_sympow_o2`);