(* We start by writing the ansatz *)
Celine[f_,n_,k_,r_,s_]:=
Module[{ansatz},
ansatz=Sum[c[i,j]*(f/.{n->n+i,k->k+j}),{i,0,r},{j,0,s}]
];
Celine[Binomial[n,k],n,k,1,1]

(* We have to divide this ansatz by f and simplify. *)
Celine[f_,n_,k_,r_,s_]:=
Module[{ansatz},
ansatz=FunctionExpand[Sum[c[i,j]*(f/.{n->n+i,k->k+j}),{i,0,r},{j,0,s}]/f]
];
Celine[Binomial[n,k],n,k,1,1]

(* Let's try a more complicated example *)
expr=Product[(Table[RandomInteger[{-10,10}],{3}].{n,k,1})!^((-1)^i),{i,6}]

(* This is quite slow *)
Timing[Celine[expr,n,k,1,1];]

(* Better perform the simplification on each part separately *)
Celine[f_,n_,k_,r_,s_]:=
Module[{ansatz},
ansatz=Sum[c[i,j]*FunctionExpand[(f/.{n->n+i,k->k+j})/f],{i,0,r},{j,0,s}]
];
Celine[Binomial[n,k],n,k,1,1]

(* Much better. *)
Timing[Celine[expr,n,k,1,1];]

(* But what about this? :-( *)
Timing[Celine[expr,n,k,10,10];]

(* Using the rational functions (certificates) of the hypergeometric input *)
Celine[f_,n_,k_,r_,s_]:=
Module[{u,v,ansatz},
u=FunctionExpand[(f/.n->n+1)/f];
v=FunctionExpand[(f/.k->k+1)/f];
ansatz=Table[Product[(u/.{n->n+i1,k->k+j}),{i1,0,i-1}]*Product[(v/.k->k+j1),{j1,0,j-1}],{i,0,r},{j,0,s}]
];

(* This is acceptable *)
Timing[Celine[expr,n,k,10,10];]

(* Basic version of Sister Celine's algorithm *)
Celine[f_,n_,k_,r_,s_]:=
Module[{u,v,mat,ns,rec},
u=FunctionExpand[(f/.n->n+1)/f];
v=FunctionExpand[(f/.k->k+1)/f];
mat=Flatten[Table[Product[(u/.{n->n+i1,k->k+j}),{i1,0,i-1}]*Product[(v/.k->k+j1),{j1,0,j-1}],{i,0,r},{j,0,s}]];
(* Multiply with the common denominator *)
mat=Together[mat*(PolynomialLCM@@Denominator[mat])];
(* Coefficient comparison w.r.t. k *)
mat=Transpose[PadRight[CoefficientList[#,k]&/@mat]];
(* Compute the kernel *)
ns=NullSpace[mat];
If[ns==={},Return[{}]];
(* Each kernel vector gives a k-free recurrence. *)
(* Sum these by replacing f[n+i,k+j] by SUM[n+i]. *)
rec=ns.Flatten[Table[SUM[n+i],{i,0,r},{j,0,s}]];
(* Simplify the obtained recurrences and return them. *)
rec=Collect[Numerator[Together[#]],SUM[_],Expand]&/@rec;
Return[rec];
];

(* It works! *)
Celine[Binomial[n,k],n,k,1,1]

Celine[Binomial[n,k]^2,n,k,2,2]

FullSimplify[%[[1]]/.SUM[n_]:>Binomial[2n,n]]

(* Sister Celine's algorithm with (a quick-and-dirty version of) Verbaeten completion *)
Celine[f_,n_,k_,r_,s_,e_]:=
Module[{u,v,mat,den,deg,idx,ns,rec},
u=FunctionExpand[(f/.n->n+1)/f];
v=FunctionExpand[(f/.k->k+1)/f];
mat=Table[Product[(u/.{n->n+i1,k->k+j}),{i1,0,i-1}]*Product[(v/.k->k+j1),{j1,0,j-1}],{i,0,r+e},{j,0,s+e}];
den=PolynomialLCM@@Denominator[Flatten[mat[[1;;r+1,1;;s+1]]]];
mat=Together[den*mat];
deg=Max[Exponent[Flatten[mat[[1;;r+1,1;;s+1]]],k]];
mat=Flatten[mat];
idx=Select[Range[Length[mat]],Denominator[mat[[#]]]===1&&Exponent[mat[[#]],k]≤deg&];
mat=mat[[idx]];
mat=Transpose[PadRight[CoefficientList[#,k]&/@mat]];
ns=NullSpace[mat];
If[ns==={},Return[{}]];
rec=ns.(Flatten[Table[SUM[n+i],{i,0,r+e},{j,0,s+e}]][[idx]]);
rec=Collect[Numerator[Together[#]],SUM[_],Expand]&/@rec;
Return[rec];
];

(* Compare: original Celine's algorithm *)
Timing[Celine[(-1)^k*Binomial[2n,n+k]^2,n,k,2,4]]

(* and with Verbaeten completion *)
Timing[Celine[(-1)^k*Binomial[2n,n+k]^2,n,k,1,3,1]]

(* Another example: Apery numbers *)
Timing[Celine[Binomial[n,k]^2*Binomial[n+k,k]^2,n,k,4,3]]
Timing[Celine[Binomial[n,k]^2*Binomial[n+k,k]^2,n,k,2,3,2]]

(* Note that we do not get the minimal second-order recurrence. *)
Arec=(1+3 n+3 (n^2)+(n^3)) SUM[n]+(-117-231 n-153 (n^2)-34 (n^3)) SUM[1+n]+(8+12 n+6 (n^2)+(n^3)) SUM[2+n]

(* However, we can check that the larger recurrence is a consequence of the minimal one: *)
Together[%%[[2,1]]//.SUM[n+i_/;i≥2]:>Solve[(Arec/.n->n+i-2)==0,SUM[n+i]][[1,1,2]]]