Conformal blocks via recursion relations (1307.6856)

ν=1/2;
C1[l_,k_]:=-((k((((2 k)!))^2))/((2^(4 k-1))(((k!))^4)))(Pochhammer[l+2ν,2k]/Pochhammer[l+ν,2k]);
C2[l_,k_]:=-((k(ν+l-k))/((ν+l+k)(((k!))^2)))Pochhammer[ν,k]Pochhammer[1-ν,k](((Pochhammer[((ν+l+1-k)/2),k])^2)/((Pochhammer[((ν+l-k)/2),k])^2));

C3[l_,k_]:=-((k((((2 k)!))^2))/((2^(4 k-1))(((k!))^4)))(Pochhammer[1+l-2k,2k]/Pochhammer[1+ν+l-2k,2k]);

P[l_,η_]:=((l!)/(Pochhammer[2ν,l]))GegenbauerC[l,ν,η] ;

(*This is a function which computes h recursively to order (r^nmax).
Notice that h_here = h_1307.6856 times , i.e. the denominator of their eq (3.12). This denominator does not partecipate in the recursion and it makes sense to factor it out. *)

h[Δ_,l_,nmax_]:=P[l,η]+
Sum[Δi=1-l-2k;li=l+2k;((C1[l,k](r^(2k)))/(Δ-Δi)) h[Δi+2k,li,nmax-2k],{k,1,nmax/2}]+
Sum[Δi=1+ν-k;li=l;((C2[l,k](r^(2k)))/(Δ-Δi)) h[Δi+2k,li,nmax-2k],{k,1,nmax/2}]+
Sum[Δi=1+l+2ν-2k;li=l-2k;((C3[l,k](r^(2k)))/(Δ-Δi)) h[Δi+2k,li,nmax-2k],{k,1,Min[nmax/2,l/2]}];

Δϵ=1.41;
Hϵ[r_,η_]:=Evaluate[h[Δϵ,0,12]//Simplify]

Hϵ[r,η]

R[x_,y_]:=(Sqrt[(x^2)+(y^2)]);
Eta[x_,y_]:=Cos[ArcTan[y/x]];

Show[GraphicsRow[{Plot3D[{Evaluate[(r^(Δϵ))/.{r->R[x,y],η-> Eta[x,y]}]},{x,-1,1},{y,-1,1},AxesLabel-> {"r","η"},RegionFunction->Function[{x,y,z},x^2+y^2≤1],PlotRange-> {0,3}],Plot3D[{Evaluate[(((r^(Δϵ))Hϵ[r,η])/((((1-(r^2)))^ν)(Sqrt[(((1+(r^2)))^2)-4(r^2)(η^2)])))/.{r->R[x,y],η-> Eta[x,y]}]},{x,-1,1},{y,-1,1},AxesLabel-> {"r","η"},RegionFunction->Function[{x,y,z},x^2+y^2≤0.9],PlotRange-> {0,3}]}],ImageSize->Full]

HT[r_,η_]:=Evaluate[h[3.,2,12]//Simplify]

HT[r,η]//Expand//Collect[#,r,Simplify]&

Show[GraphicsRow[{Plot3D[{Evaluate[(r^3)P[2,η]/.{r->R[x,y],η-> Eta[x,y]}]},{x,-1,1},{y,-1,1},AxesLabel-> {"r","η"},RegionFunction->Function[{x,y,z},x^2+y^2≤1],PlotRange-> {-1,1.5}],Plot3D[{Evaluate[(((r^3)HT[r,η])/((((1-(r^2)))^ν)(Sqrt[(((1+(r^2)))^2)-4(r^2)(η^2)])))/.{r->R[x,y],η-> Eta[x,y]}]},{x,-1,1},{y,-1,1},AxesLabel-> {"r","η"},RegionFunction->Function[{x,y,z},x^2+y^2≤0.9],PlotRange-> {-1,1.5}]}],ImageSize->Full]