clear;
// lambda= input("Entre com o modulo do coeficiente de decaimento:") // 0.001
lambda = 0.001
N      = 1/lambda; // # de particula total
r_max  = 100; // # max de realizacoes
t_max  = 300; // # tempo maximo
clf();
rand('uniform');
rand('seed',2)
t1 = 0:t_max;
for r = 1:r_max
 n         = N;      // # de particulas 
 nd        = 0.;     // # de particulas decaidas
 nd_t(1,r) = [nd];   // inicializa curva
 for t = 1:t_max     // laco no tempo
   if rand() < (n/N) then // n/N propocao de part. nao decaidas
     n  = n - 1;           // decaimento de uma particula
     nd = nd + 1;          // praticula decaida
   end;
   nd_t(t+1,r) = -nd/N; // constroe a curva n_d/N versus tempo,para o decaimento 1/N <0
 end
 scf(0);
 plot2d(t1,[nd_t(:,r)],style=[-9]); // grafica as curvas
end
m_nr_t  = zeros(1,t_max+1); // curva de valor medio
m_nr2_t = zeros(1,t_max+1); // curva de valor ao quadrado medio
m_sup   = zeros(1,t_max+1); // curva de valor ao quadrado medio
m_inf   = zeros(1,t_max+1); // curva de valor ao quadrado medio
for  t   = 1:t_max+1
 for r   = 1:r_max
   m_nr_t(t)  = m_nr_t(t) + nd_t(t,r);
   m_nr2_t(t) = m_nr2_t(t) + nd_t(t,r)*nd_t(t,r);
 end
 m_nr_t(t)  = m_nr_t(t)/(r_max);                            // valor medio
 m_nr2_t(t) = m_nr2_t(t)/(r_max);                           // valor quadrado medio
 m_nr2_t(t) = sqrt(abs(m_nr2_t(t) - m_nr_t(t)*m_nr_t(t)));  // desvio padrao
 m_inf(t)   = m_nr_t(t)           - m_nr2_t(t)/sqrt(r_max); // media - erro-padrao
 m_sup(t)   = m_nr_t(t)           + m_nr2_t(t)/sqrt(r_max); // media + erro-padrao
end                                                         // erro-padrao = desvio-padrao/raiz(realizacoes)
//clf();
scf(1);
plot([t1],[m_nr_t],[t1],[m_inf],[t1],[m_sup]);//figure(2)
legend("valor medio","media - erro-padrao","media + erro-padrao");//legenda