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

(*The lattice defined in (2.1)*)
xx[s_]:=s^2+beta1 s
yy[t_]:=t^2+beta2 t

(*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) (-1+k-n) (beta0 (6+k+3 beta1 k-2 (k^2)-13 n-3 beta1 n+4 n (k+n)-6 (Nn^2)-6 beta3 (1-n+Nn))+(k-n) (beta3 (-7+2 k+4 n-6 Nn)+2 ((-1+n) (-4+2 k+n)-3 (Nn^2)))+beta1 (beta3 (-3 (-2+k+n)+6 Nn)+6 (-1+k+n-k n+(Nn^2))))

s2[n_,k_]:=-(1/6) k (beta3 (5+(3-2 k) k+6 (-2+n) n+3 beta2 (-1+k-2 Nn)+6 (1+k-2 n) Nn)-2 ((-1+n) (-2-3 beta2 k+2 (k^2)+3 (1+beta2-k) n)+3 (-1+beta2-k+2 n) (Nn^2))+beta0 ((-1+k) (5+2 k-6 n)-3 beta2 (1+k-2 n)+6 (Nn^2)+6 beta3 (1-n+Nn)))

(*The coefficients of the matrix T_n*)
t1[n_,k_]:=-(1/180) (-1+k-n) (k-n) (-15 (beta1^2) (-12-3 (-3+k) k+7 n-2 k n-(n^2)+6 (Nn^2)+beta3 (8-4 k-2 n+6 Nn))+30 beta1 (-(-1+k) (-1+k (-5+2 k)-2 (-3+n) n)+(1+4 k-4 n) (Nn^2)+beta3 (1-3 (k^2)+(-5+n) n+2 k (2+n)+(1+4 k-4 n) Nn))+6 beta0 (4+34 k-26 (k^2)+4 (k^3)+5 (beta1^2) (1+k-n)+(29-8 k) (n^2)-4 (n^3)-5 (Nn^2)-20 k (Nn^2)-5 beta3 (1+(5-2 k) k+2 (-3+n) n+(1+4 k-4 n) Nn)+2 n (-22+k+4 (k^2)+10 (Nn^2))+5 beta1 (-3-2 (k^2)+k (7-2 n)+(-1+n) n+3 (Nn^2)+3 beta3 (1-k+Nn)))+(1+k-n) (-18+145 k-133 n+2 (2 (k-n) (5 (k^2)+2 (-9+n) n+k (-27+8 n)-15 (Nn^2))+3 beta3 (1+6 (k^2)+n (15-4 n+10 Nn)-k (15+2 n+10 Nn)))))

t2[n_,k_]:=(1/18) k (k-n) (3 beta1 (-1+5 k+2 n+beta3 (2-2 k (3+k-3 n)-3 n-6 (-1+n) Nn)-2 (k (k+(k^2)+3 n-2 k n)+3 (-1+n) (Nn^2))-3 beta2 (-1+k (-k+n)+2 (Nn^2)+beta3 (1-k+2 Nn)))+(1+k-n) (19-24 beta3-19 k+9 beta3 k-4 (k^2)+6 beta3 (k^2)+4 (k^3)-6 (-1+beta3+k) (n^2)-30 beta3 Nn-6 beta3 k Nn-30 (Nn^2)-6 k (Nn^2)+2 n (-13+15 k-2 (k^2)+12 (Nn^2)+3 beta3 (5-2 k+4 Nn))+3 beta2 (1-2 (-2+k) k+2 (-3+n) n+6 (Nn^2)+beta3 (3-3 k+6 Nn)))+3 beta0 (8-10 beta3-5 k-5 (k^2)+4 beta3 (k^2)+2 (k^3)+3 beta1 (1+beta2-beta3-k) (1+k-n)-4 (-1+beta3+k) (n^2)-12 beta3 Nn-6 beta3 k Nn-12 (Nn^2)-6 k (Nn^2)+beta2 (2+6 k-9 n-2 (k-n) (k+2 n)+6 (Nn^2)+beta3 (3-3 k+6 Nn))+n (-13+13 k+12 (Nn^2)+beta3 (15-4 k+12 Nn))))

t3[n_,k_]:=-(1/180) k (1+k) (-30 (n^2) (beta3+2 k-4 beta3 k-2 (k^2)+6 beta3 Nn+6 (Nn^2))-15 (beta2^2) (-3 (4+k (3+k))+8 (2+k) n-6 (n^2)+6 (Nn^2)+beta3 (2-2 k+6 Nn))+(2+k) ((-1+k) (9+4 k (17+5 k))-120 (Nn^2)-12 beta3 (1-k (9+2 k)+10 Nn))+12 n (1+14 k-3 (k^2) (3+2 k)+5 (7+2 k) (Nn^2)+5 beta3 (1-k (5+2 k)+(7+2 k) Nn))+30 beta2 (-(1+k-n) (-1+k (5+2 k-4 n)+n)+(7+2 k-6 n) (Nn^2)+beta3 (-2 (k^2)-3 (((-1+n))^2)+k (-7+6 n)+(7+2 k-6 n) Nn))+6 beta0 (-1-14 k+9 (k^2)+6 (k^3)+5 (beta2^2) (4+2 k-3 n)-35 (Nn^2)-10 k (Nn^2)-5 beta3 (1-n+k (-5-2 k+4 n)+(7+2 k-6 n) Nn)+10 n (k-(k^2)+3 (Nn^2))+5 beta2 (1-n+k (-4-3 k+4 n)+3 (Nn^2)+3 beta3 (-1-k+n+Nn))))

(*We first check the recurrence of the first family of bivariate Racah defined in (1.6)*)

(*Univariate Racah defined in (1.1)*)
rr[alpha_,beta_,gamma_,delta_,s_,n_]:=Pochhammer[alpha+1,n]Pochhammer[beta+delta+1,n]Pochhammer[gamma+1,n]HypergeometricPFQ[{-n,n+alpha+beta+1,-s,s+gamma+delta+1},{alpha+1,beta+delta+1,gamma+1},1]

(*First family of Bivariate Racah defined in (1.6)*)
RR1[s_,t_,beta0_,beta1_,beta2_,beta3_,Nn_,n_,m_]:=rr[beta1-beta0-1,beta2-beta1-1,-t-1,beta1+t,s,n]rr[2n+beta2-beta0-1,beta3-beta2-1,n-Nn-1,n+beta2+Nn,t-n,m]

(*The polynomial vector of polynomials RR1[...,n-k,k]*)
PP1[n_]:=Table[RR1[s,t,beta0,beta1,beta2,beta3,Nn,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 Pochhammer[beta1-beta0+n-l,l-k]Pochhammer[-k,k]Pochhammer[k-n,n-l]Pochhammer[beta2-beta0-1+n-k,n-l]Pochhammer[2n-k+beta3-beta0-1,k])/((n-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(beta3-beta0+n-1)

(*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 for the first family*)
vvv1[n_]:=Factor[xx[s]PP1[n]-AAA1[n,1].PP1[n+1]-BBB1[n,1].PP1[n]-CCC1[n,1].PP1[n-1]]

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

vvv1[1]

vvv1[2]

vvv1[3]

vvv2[1]

vvv2[2]

vvv2[3]

(*We check the recurrence of the second family of bivariate Racah defined in (2.8)*)

(*Univariate Racah defined in (1.1)*)
rr[alpha_,beta_,gamma_,delta_,s_,n_]:=Pochhammer[alpha+1,n]Pochhammer[beta+delta+1,n]Pochhammer[gamma+1,n]HypergeometricPFQ[{-n,n+alpha+beta+1,-s,s+gamma+delta+1},{alpha+1,beta+delta+1,gamma+1},1]

(*Second family of Bivariate Racah defined in (2.8)*)
RR2[s_,t_,beta0_,beta1_,beta2_,beta3_,Nn_,n_,m_]:=rr[2m-beta1+beta3-1,beta1-beta0-1,m-Nn-1,m-Nn-beta1,Nn-m-s,n]rr[beta3-beta2-1,beta2-beta1-1,s-Nn-1,-beta2-Nn-s,Nn-t,m]

(*The polynomial vector of polynomials RR2[...,n-k,k]*)
PP2[n_]:=Table[RR2[s,t,beta0,beta1,beta2,beta3,Nn,n-k,k],{k,0,n}]

(*The matrix G_n defined in (3.8)*)
G2[n_]:=Table[g2[n,k,l],{k,0,n},{l,0,n}]

g2[n_,k_,l_]:=Binomial[k,l]Pochhammer[beta2-beta3-k+1,k-l]Pochhammer[k-beta1+beta3-1,l]Pochhammer[k+n-beta0+beta3-1,n-k]

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

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

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

lambda[n_]:=n(beta3-beta0+n-1)

(*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*)
AAA2[n_,j_]:=G2[n].LL[n,j].Inverse[G2[n+1]]
BBB2[0,j_]:=(-AAA2[0,j].GG2[1]).Inverse[G2[0]]
BBB2[n_,j_]:=(GG2[n].LL[n-1,j]-AAA2[n,j].GG2[n+1]).Inverse[G2[n]]
CCC2[1,j_]:=(-AAA2[1,j].GGG2[2]-BBB2[1,j].GG2[1]).Inverse[G2[0]]
CCC2[n_,j_]:=(GGG2[n].LL[n-2,j]-AAA2[n,j].GGG2[n+1]-BBB2[n,j].GG2[n]).Inverse[G2[n-1]]

(*The recurrence for the second family*)
uuu1[n_]:=Factor[xx[s]PP2[n]-AAA2[n,1].PP2[n+1]-BBB2[n,1].PP2[n]-CCC2[n,1].PP2[n-1]]

uuu2[n_]:=Factor[yy[t]PP2[n]-AAA2[n,2].PP2[n+1]-BBB2[n,2].PP2[n]-CCC2[n,2].PP2[n-1]]

uuu1[1]

uuu1[2]

uuu1[3]

uuu2[1]

uuu2[2]

uuu2[3]