esp=0;
Esp2=0;
var=0;
n=10000;
lambda=4;


for i = 1:n;
    u=rand(i);
    cumul=exp(-lambda);
    ajout=cumul;
    alpha=0;
    while (cumul<u) do
        alpha = alpha + 1;
        ajout = (ajout * lambda)/alpha
        cumul = cumul + ajout;
    else
        poisson = alpha;
    end

    esp = esp + (poisson/n);
    Esp2 = (poisson^2)/n + Esp2;
end;
var = Esp2 - esp^2;

disp(var);
disp(esp);
