CalcolaB=function(P){
     X = diag(nrow(P))
     Y = X
     massima = nrow(P)-1
     for(potenza in 1:massima)
     {
         X <- X %*% P
         Y <- Y + X
         return(Y)
     }
 }


====================================================


Markov=function(P, nrepl, passi){
    prova =1 
    Y = vector( length = passi)
    for(prova in 1:nrepl)
    {
        xi=runif(passi, min=0, max=1)
        dim(xi)=c(passi,1)
        passo = 0
        i <- sample(1:nrow(P), size = 1)
        print(paste("Sorteggio X[0](omega)", i))
        X = matrix(, ncol = passi, nrow = nrepl)
        arresto = 0
        k=1
        j=1
        while(arresto < xi[1])
        {
            arresto = arresto + P[i,j]
            X[prova,1] = j
            Y[1] = X[prova,1]
            j= j+1
        }
                for(k in 2:passi)
        {
            i=X[prova,k-1]
            arresto = 0 
            j=1
            while(arresto < xi[k])
            {
                arresto = arresto + P[i,j]
                X[prova,k] = j
                j=j+1
            }
            Y[k] <- X[prova,k]
        }
        print(Y)   
        rosso = prova/nrepl
        blu = 1 - rosso
        plot(Y, col = rgb(rosso,0,blu))
    }
}