esp=0;
Esp2=0;
var=0;
alph=0;
n=10000;
lambda=4;
ajout=1;


for i = 1:n;
    somme=0;
    u=rand(i);
        if (u>somme) then
            alph = alph + 1;
            ajout = (ajout * lambda)/alph;
            somme = somme + ajout;
        end;
    esp = esp + (alph/n);
    Esp2 = (alph^2)/n + Esp2;
end;
var = Esp2 - esp^2;
esp = esp/1200;
var = var/2008839;
disp(var);
disp(esp);
