function [S N] = poisson(lambda, Tmax)

S(1) = exprnd(lambda);  
i = 1;
% continuo ad estrarre nuovi incrementi finché non ho superato la lunghezza
% dell'intervallo
N=0;    j=0;
while S(i)<Tmax
    i = i+1;
    S(i) = S(i-1) + exprnd(lambda);
    j = j + 1;
    N=[N, j]; % aggiungi evento alla lista N
end

S(end) = Tmax;
N(end) = N(end-1);