H := binomial(3*n,3*k)^2*binomial(3*n,3*k+1);

R1 := simplify(convert(subs(k = k+1, H)/H, GAMMA));

dn := (3*k+1)*(3*k-1)^3*k^3*(3*k-2)^2;

in_Delta_k := collect(normal( subs(k=k+1, dn) * R1) - dn, k, factor);

phi_subs :=  k = n-k-1;
	# The -1 is to simplify the ad-hoc reduction, shouldn't be
	# needed in a more systematic implementation.

phi_H := subs(phi_subs, add(u[i]*k^i, i=0..8)) * simplify(convert(subs(phi_subs, H)/H,GAMMA));
Ansatz := add(c[i]*k^i, i=0..8);
normal(phi_H - subs(k = k+1, Ansatz) * R1):
factor(solve({coeffs(rem(numer(%), denom(%), k), k)}, {seq(c[i],i=0..8)})):
subs(%, Ansatz):
subs(k=k+1, %) * R1 - %:
phi_H := collect(normal(phi_H - %), k, normal);
degree(%,k);

M := [seq(collect(coeff(phi_H, u[j]), k, normal), j=0..8)];
M := Matrix([seq([seq(coeff(v,k,j),j=0..8)], v=M)]);

BasisNplus  := [seq( (k^i + coeff(phi_H, u[i]))/2, i = 0..4)];
BasisNminus := [seq( (k^i - coeff(phi_H, u[i]))/2, i = 0..3)];

phi3 := add(u[i]*(k + 1/3)^i, i=0..8) * simplify( convert(subs(k=k+1/3, H)/H, GAMMA) );
Ansatz := add(c[i]*k^i, i=0..8);
normal(phi3 - subs(k = k+1, Ansatz) * R1):
factor(solve({coeffs(rem(numer(%), denom(%), k), k)}, {seq(c[i],i=0..8)})):
subs(%, Ansatz):
subs(k=k+1, %) * R1 - %:
phi3 := collect(normal(phi3 - %), k, normal);
degree(%,k);

M3 := [seq(collect(coeff(phi3, u[j]), k, normal), j=0..8)];
M3 := Matrix([seq([seq(coeff(v,k,j),j=0..8)], v=M3)]);

map(normal, M^2);
map(normal, M3^3);

# Here are bases of the 4 submodules:

Npp := LinearAlgebra:-NullSpace( Matrix( [[ M3 - 1 ], [ M-1 ]]) );
Npm := LinearAlgebra:-NullSpace( Matrix( [[ M3 - 1 ], [ M+1 ]]) );
Nmp := LinearAlgebra:-NullSpace( Matrix( [[ M3^2+M3+1 ], [ M-1 ]]) );
Nmm := LinearAlgebra:-NullSpace( Matrix( [[ M3^2+M3+1 ], [ M+1 ]]) );

# The dimensions of these submodules match the orders of the operators Lpp, Lmp, Lpm, Lmm
# computed in the other file:

nops(Npp);
nops(Nmp);
nops(Npm);
nops(Nmm);

# We can also compute Lpp, Lmp, Lpm, Lmm by projecting H onto each of these 4 subspaces
# and computing a telescoper for each projection, as in the binomial(n,k)^/(2*n+3*k) example.