Mdn[x_, rho_, n_] :=
Module[{i = 0, I = 0, GuessP, Sum1, P, S, Difference, Max = 0},
  
  GuessP[0] = ( Sum[(rho n)^i/i!,{i,0,n-1}] + (rho n)^n/n! 1/(1-(rho n)/n) )^-1;
  For[i=1, i 10^-12),
      i=i+1;
      GuessP[i] = GuessP[i-1] (rho n)/n;
  ];
  I = i; 
  For[i=(I+1), i<=(I+n), i++,
      GuessP[i] = GuessP[i-1] (rho n)/n;
  ];

  While[ ((Max > 10^-12) || (Max == 0)),
      Max=0;
      Sum1 = Sum[GuessP[k], {k,0,n}];
      For[i=0, i<=I, i++,
	  P[i] = Sum1 (rho n)^i/i! E^-(rho n)
	  + Sum[GuessP[k] (rho n)^(n+i-k)/(n+i-k)! E^-(rho n), {k,n+1,n+i}];
      ];
      For[i=(I+1), i<=(I+n), i++,
	  P[i] = P[i-1] (rho n)/n;  
      ];
      S = Sum[P[k], {k,0,I+n}];
      For[i=0, i<=(I+n), i++,
	  P[i] = P[i]/S;
      ];
      For[i=0, i<=(I+n), i++,
	  Difference = Abs[P[i] - GuessP[i]];
	  If [Difference > Max,
	      Max = Difference
	  ];
	  GuessP[i] = P[i]; 
      ];
  ];
  P[i_]:=0;
  Table[P[i],{i,0,x}]
]