Clear[xx,yy,s,t,s1,s2,t1,t2,t3,n,k,beta0,beta1,beta2,beta3,Nn,S,T,PP,g1,G1,GG1,GGG1,AAA1,BBB1,CCC1,g2,G2,GG2,GGG2,AAA2,BBB2,CCC2,Z,lambda,l,LL,i,j,vvv1,vvv2,uuu1,uuu2,x,y]

(*The lattice*)
xx[s_]:=s^2
yy[t_]:=t^2

(*The matrix S_n defined in (3.5)*)
S[n_]:=Table[If[k==l,s1[n,k],If[k==l+1,s2[n,l],0]],{k,1,n+1},{l,1,n}]

(*The matrix T_n defined in (3.6)*)
T[n_]:=Table[If[k==l,t1[n,l],If[k==l+1,t2[n,l],If[k==l+2,t3[n,l],0]]],{k,1,n+1},{l,1,n-1}] 

(*The coefficients of the matrix S_n*)
s1[n_,k_]:=(1/6)(n-k+1)(4n^2+4k n-13n-2k^2-5k+12)+(n-k+1)(n-k)(a+b+c+2e2)+(k-1)(n-k+1)(1+2a+b+c+2e2)+(n-k+1)(a b+a c+b c+(2a+b+c)e2+e2^2)

s2[n_,k_]:=-(1/6)(k-1)k(-6n+2k+5)+k(n-k)(b+c)+k(k-1)(a+b+c+e2)+k(a b+a c+b c+(b+c)e2)

(*The coefficients of the matrix T_n*)
t1[n_,k_]:=-(1/30)(n-k)(n-k+1)(-9-49 k+46 (k^2)-4 (k^3)+64 n-22 k n-8 (k^2) n-29 (n^2)+8 k (n^2)+4 (n^3))-(n-k)(n-k+1)a(b+e2)(c+e2)-(k-1)(n-k)(n-k+1)a(1+b+c+2e2)-(1/3)(n-k-1)(n-k)^2(n-k+1)(a+b+c+2e2)-((k-1)/6)(n-k)(n-k+1)(4n-4k-1)(1+2a+b+c+2e2)-(1/6)(n-k)(n-k+1)(4n-4k-1)(a b+a c+b c+(2a+b+c)e2+e2^2)-(k-2)(k-1)(n-k)(n-k+1)a

t2[n_,k_]:=-(1/6)(k-1)k(n-k)(4n^2-19n+19+5 k-2 (k^2))-(k/6)(n-k)(n-k-1)(4n-4k-5)(b+c)-k(k-1)(n-k)(n-k-1)(1+a+2b+2c+2e2)-(1/6) (-1+k) k (-5+4 k)(n-k)(1+2a+b+c+2e2)-k(n-k)(n-k-1)(b+a b+c+a c+2b c+2b e2+2c e2)+k(k-1)(n-k)(-a-b-2 a b-c-2 a c-b c-e2-2 a e2-2 b e2-2 c e2-(e2^2))+k(n-k)(-a b-a c-b c-2 a b c+(-b-2 a b-c-2 a c-2 b c) e2+(-b-c) (e2^2))

t3[n_,k_]:=(1/30) (-1+k) k (1+k) (1+15 k+6 (k^2)-10 k n)-k(k+1)(n-k-1)b c-(1/6) k (1+k) (-1+4 k)(n-k-1)(b+c)-k(k+1)b c(a+e2)-(1/6) k (1+k) (-1+4 k)(a b+a c+b c+(b+c)e2)-(1/3) (-1+k) (k^2) (1+k)(a+b+c+e2)

(*Univariate Continuous dual Hahn*)
cdualHahn[a_,b_,c_,x_,n_]:=Pochhammer[a+b,n]Pochhammer[a+c,n]HypergeometricPFQ[{-n,a+I x,a-I x},{a+b,a+c},1]

(*Bivariate continuous dual Hahn defined in (4.10)*)
bicdualHahn[a_,b_,c_,e2_,x_,y_,n_,m_]:=cdualHahn[a,e2+I y,e2-I y,x,n]cdualHahn[n+a+e2,b,c,y,m]

(*The polynomial vector of polynomials bicdualHahn[...,n-k,k]*)
PP[n_]:=Table[bicdualHahn[a,b,c,e2,x,y,n-k,k],{k,0,n}]

(*The matrix G_n defined in (3.7)*)
G1[n_]:=Table[g1[n,k,l],{k,0,n},{l,0,n}]

g1[n_,k_,l_]:=(-1)^(n-k-l)Binomial[n-k,l-k]

(*The matrix G_{n,n-1} defined in (3.4)*)
GG1[n_]:=G1[n].S[n].Inverse[Z[n-1,n]]

(*The matrix G_{n,n-2} defined in (3.4)*)
GGG1[n_]:=(G1[n].T[n]+GG1[n].S[n-1]).Inverse[Z[n-2,n]]

Z[n_,l_]:=(lambda[n]-lambda[l])IdentityMatrix[n+1]

lambda[n_]:=n

(*The matrices defined in (3.3)*)
LL[n_,1]:=If[n≤-1,0,Table[If[k==l,1,0],{k,1,n+1},{l,1,n+2}]]
LL[n_,2]:=If[n≤-1,0,Table[If[k+1==l,1,0],{k,1,n+1},{l,1,n+2}]]

(*The coefficients appearing in Theorem 8*)
AAA1[n_,j_]:=G1[n].LL[n,j].Inverse[G1[n+1]]
BBB1[0,j_]:=(-AAA1[0,j].GG1[1]).Inverse[G1[0]]
BBB1[n_,j_]:=(GG1[n].LL[n-1,j]-AAA1[n,j].GG1[n+1]).Inverse[G1[n]]
CCC1[1,j_]:=(-AAA1[1,j].GGG1[2]-BBB1[1,j].GG1[1]).Inverse[G1[0]]
CCC1[n_,j_]:=(GGG1[n].LL[n-2,j]-AAA1[n,j].GGG1[n+1]-BBB1[n,j].GG1[n]).Inverse[G1[n-1]]

(*The recurrence*)
vvv1[n_]:=Factor[xx[x]PP[n]-AAA1[n,1].PP[n+1]-BBB1[n,1].PP[n]-CCC1[n,1].PP[n-1]]

vvv2[n_]:=Factor[yy[y]PP[n]-AAA1[n,2].PP[n+1]-BBB1[n,2].PP[n]-CCC1[n,2].PP[n-1]]

vvv1[1]

vvv1[2]

vvv1[3]

vvv1[4]

vvv2[1]

vvv2[2]

vvv2[3]

vvv2[4]