(* the iterations we want to perform *)
f[M_]:=Simplify[Inverse[ Inverse[{{1,T},{0,1} }. M . {{1,0},{T,1}}] + (1/r){{1,0},{0,0}}]]

(* iterate the first time *)
P1=f[(1/eps){{1,0},{0,1}}];
MatrixForm[P1]

P1 ( 1 + eps r + (T^2))//Simplify//MatrixForm

(* iterate a second time *)
P2=f[f[(1/eps){{1,0},{0,1}}]]

(* lets take eps equal to zero now *)
P2p = P2 /.{eps->0};
MatrixForm[P2p]

(* iterate a third time *)
f[P2p]//MatrixForm

(* a fourth time ... we begin to see the pattern *)
f[f[P2p]]//MatrixForm

(* Lets do the general case *)
Clear[Pkp1];
Pkp1[k_]:=2 (r/((k+1)(k+2))){{2k+1,(3/T)},{(3/T),(6/(k (T^2)))}};
MatrixForm[Pkp1[k]]

value=f[Pkp1[k]]

(* Does this equal Pkp1[k+1]? Yes*)
value - Pkp1[k+1]//Simplify