restart;
 N:=5:
 a[0,0]:=1:psi[0]:=2*Y-2:
 for p from 1 to N do a[p,p]:=0 end do:
  for p from 0 to N do
    for t from 0 to N do
     if (p=0) and (t=0) then else
       K[1]:=simplify(add(add(b[k,i]*(X+p-t)/(t-i)
        *expand(binomial(X+p-i-1,p-k-1))*binomial(p-k-1,p+i-k-t)
        *(-1)^(p+k)*2^i,i=max(0,k+t-p)..t-1),k=0..p-1)):
     K[2]:=simplify(add(b[k,t]
        *expand(binomial(X+p-t,p-k))*(-1)^(p+k)*2^t,k=0..p-1)):
     K[3]:=simplify(add(add(add(chi[j]*b[k,i]
        *expand(binomial(X+p-i-j,p-j-k))*expand(binomial(X+k-i-1,t-i-j))
        *(-1)^(k+i)*2^(i+j-1),i=0..t-j),j=0..min(p-k,t)),k=0..p)):
    if t<p then
      a[p,t]:=simplify(solve(K[1]+K[2]+K[3]=(2^p*Y-2^t)*b[p,t],b[p,t])):
    elif t=p then
      psi[p]:=simplify(solve(K[1]+K[2]+K[3]=0,chi[p])):
    elif t>p then
     a[p,t]:=simplify(solve(K[1]+K[2]+K[3]=(2^p*Y-2^t)*b[p,t],b[p,t]))
    end if:
    end if:
   chi[p]:=psi[p]:
   b[p,t]:=a[p,t]:
 end do:
 end do:
 psi[N];