RecMd1[x_, rho_] :=
Module[{n=1, i=1, P},
  P[0] = 1-rho;
  P[1] = (1-rho) (E^rho-1);
  For[i=1, i < x, i++,
      P[i+1] = 1/(E^-rho) ( P[i] - (P[0]+P[1]) rho^i/i! E^-rho
               - Sum[P[n] rho^(i-n+1)/(i-n+1)! E^-rho,{n,2,i}] )];
  Table[P[i], {i,0,x}]
]