# WARNING: This doesn't attempt to reduce so use gfun[rec*rec] if both recs are >= 2nd order 
# 1.3 v local, no globals, and u,n are passed in
# 1.2 u,v,w,h are global so they can be changed outside of the procedure
# 1.1 returns the recurrence in terms of u (not h)
multrecs:=proc(L1, L2, u, n)  # input 2 lists of coeffs (decreasing) order, each minus u(n+ord): u(n+1)-n*u(n) inputs as [-n]
local m1, m2, max12, prod12, L1t, L2t, r1, r2, v, w, i, wcoeffs, coeffsol, c, uvs, x, eqn, eqns, weqn, h,j;
m1,m2 := nops(L1), nops(L2);
max12 := max(m1,m2);
prod12 := m1*m2;
L1t, L2t := L1[1..max(nops(L1)-1, 1)], L2[1..max(nops(L2)-1, 1)];  # Seeing if we have zero coeffs for all but first (first not in list) and last
if [L1t,L2t] = [[seq(0, i=1..nops(L1t))],[seq(0, i=1..nops(L2t))]] then prod12 := lcm(m1,m2) fi;
r1 := u(n + m1) = -add(L1[m1-i]*u(n+i), i=0..m1-1);
r2 := v(n + m2) = -add(L2[m2-i]*v(n+i), i=0..m2-1);
r1 := subs(n=n + max12-m1, r1);
r2 := subs(n=n + max12-m2, r2);
#r1 := u(n + max12) = -add(L1[max12-i]*u(n+i), i=max12-m1..max12-1);
#r2 := v(n + max12) = -add(L2[max12-i]*v(n+i), i=max12-m2..max12-1);
w[n + max12] := [rhs(r1) , rhs(r2)];
for i from max12+1 to (max12 + prod12) do
 w[n+i] := [subs(r1, subs(n=n+1, w[n+i-1][1])) , subs(r2, subs(n=n+1, w[n+i-1][2]))];
od:
wcoeffs := {seq(c[n+i], i=max12..max12+prod12-1)};
uvs := {seq(u(n+i)=x^((i+1)*prod12), i=max12-m1..max12-1), seq(v(n+i)=x^(i), i=max12-m2..max12-1)};
eqn := convert(w[n+max12 + prod12], `*`) + add(c[n+j] * convert(w[n+j], `*`), j=max12..(max12 + prod12 - 1));
for i to nops(uvs) do
 eqn := subs(uvs[i], eqn)
od;
eqn := collect(eqn, x);
eqns := { coeffs(eqn, x) };
coeffsol := solve(eqns, wcoeffs);
if coeffsol = NULL then return('Fail') fi;
#unassign(w);
weqn := h(n+max12 + prod12)+subs(coeffsol, add(c[n+j] * h(n+j), j=max12..(max12 + prod12 - 1)));
#subs(N=n, %);
subs(n=n-max12, weqn); #Just moving subscripts back down (so lowest subscript is n)
collect(%, {h(n+2), h(n+1), h(n)}, factor);
subs(h=u, %%);
end: