%/* Section 1. Preliminary definitions to facilitate studying the running example */;

%/* computations E_l, where UU and VV stand for formal parameters representing u- and v- tuples of uknown lengths;

E(l, VV, UU) := block ( [],
if (l=0) then 
    1
    else ( 
        -(1/l)*sum(E(l-k,VV, UU)*sum(binomial(k,ll)*UU(ll)*VV(k-ll), ll, 0, k), k, 1, l)
    )
);

%/* assigning actual u- and v-parameters */;

U(l):=u[l];V(l):=v[l];

%/* demo calculations */;

E(0,V,U);E(1,V,U);

expand(ratsimp(E(2,V,U)));

expand(ratsimp(E(3,V,U)));

%/* computing power sums */;

expon(x,y) := block([], 
if (y=0) then 1 else x^y
);

expon(x,0);

expon(0,0);

expon(x,1);

p(t1,t2,l):=expon(t1,l)+expon(t2,l);

/* demo calculations */;

p(1,1,0);p(1,2,0);p(2,2,0);p(2,4,0);

p(1,1,1);p(1,2,1);p(2,2,1);p(2,4,1);

p(1,1,2);p(1,2,2);p(2,2,2);p(2,4,2);

p(1,1,3);p(1,2,3);p(2,2,3);p(2,4,3);

;

%/* substituion v-parameter by the power sum */;

Vt(l):=p(t1,t2,l);

Vt(2);Vt(3);

;

%/* claculation E_l with v-parameters substituted by power-sums;

EE(l, t1, t2) := block ([],
    retVal : E(l, V, U),
    for ll:0 while ll<=l do (
      retVal:ratsubst(p(t1,t2,ll), v[ll], retVal )
    ),
    retVal
);

;

%/* demo calculations */;

EE(1, 1, 2);

EE(2, 1, 2);

%/* Encoding S_l as function on alpha-parameters (here a-parameters) */;

S(l,a) := block ( [],
    retVal:0,
for i:1 while i<=4 do (
    for j:i while j<=4 do
        (
         retVal:retVal+a(i,j)*EE(l, i, j)
        )
    ),
   retVal 
);

S(0, a);

S(1, a);

S(2, a);

%/* Section 2. Studying the running example*/;

lhs(x,P):=P(x-1)*P(x-1)-3*P(x-1)*P(x-2)+(5/2)*P(x-2)*P(x-2)-(1/2)*P(x-2)*P(x-4)-P(x)+2*P(x-1)-(1/8)*P(x-2);

;

%/* encoding alpha-s for the running example */;

ab(i,j):= block([],
if (i=1 and j=1) then 1 else
(if (i=1 and j=2) then -3 else (
        if (i=2 and j=2) then 5/2 else (
            if (i=2 and j=4) then -1/2 else 0
        ) 
        )
    )
);

%/* computing S_0, S_1, S_2 for the runnig example */;

S(0, ab); ratsimp(S(1, ab)); ratsimp(S(2, ab));

%/* applying unknow coefficients to obtain solutions of degree 3 */;

Pinst(x):=a3*x^3+a2*x^2+a1*x+a0;

lhsInst:expand(lhs(x,Pinst));

%/*creating and solving a system w.r.t. the unknown coefficients a0, a1, a2, a3*/;

C6(a3, a2, a1, a0) := ratcoeff(lhsInst, x^6); C6(a3, a2, a1, a0);

C5(a3, a2, a1, a0) := ratcoeff(lhsInst, x^5); C5(a3, a2, a1, a0);

C4(a3, a2, a1, a0) := ratcoeff(lhsInst, x^4); C4(a3, a2, a1, a0);

C3(a3, a2, a1, a0) := ratcoeff(lhsInst, x^3); C3(a3, a2, a1, a0);

C2(a3, a2, a1, a0) := ratcoeff(lhsInst, x^2); C2(a3, a2, a1, a0);

C1(a3, a2, a1, a0) := ratcoeff(lhsInst, x^1); C1(a3, a2, a1, a0);

C0(a3, a2, a1, a0) := ratcoeff(lhsInst, x,0); C0(a3, a2, a1, a0);

algsys([C3(a3,a2,a1,a0)=0, C2(a3,a2,a1,a0)=0, C1(a3,a2,a1,a0)=0, C0(a3,a2,a1,a0)=0 ], [a3, a2, a1, a0]);

%/* studying a shifted equation */;

algsys([(2+Delta)*(2+Delta) = (4+Delta)*(1+Delta)], [Delta]);

Sshift(l,a, Delta) := block ( [],
    retVal:0,
for i:1 while i<=4 do (
    for j:i while j<=4 do
        (
         retVal:retVal+a(i,j)*EE(l, i+Delta, j+Delta)
        )
    ),
   retVal 
);

lhsShifted(x,P):=P(x-1-Delta)*P(x-1-Delta)-3*P(x-1-Delta)*P(x-2-Delta)+(5/2)*P(x-2-Delta)*P(x-2-Delta)-(1/2)*P(x-2-Delta)*P(x-4-Delta)-P(x)+2*P(x-1-Delta)-(1/8)*P(x-2-Delta);

Sshift(0,ab, Delta);

ratsimp(Sshift(1,ab, Delta));

ratsimp(Sshift(2,ab, Delta));