function traj = saut_incluse(L,x0,n)
//     Simule par la chaine incluse, n sauts du processus 
//     de generateur L a partir de x0. 
//     Retourne une matrice a n+1 colonnes et 2 lignes. 
//     La premiere ligne contient les etats visites, la
//     seconde contient les durees de sejour.

ts = -diag(L);                   // taux de sejour 
P = (L+diag(ts))./(ts*ones(ts'));// matrice de transition
etats = grand(n,"markov",P,x0);  // etats visites
etats = [x0,etats];
durees = grand(etats,"exp",1);   // exponentielles
durees = durees./(ts(etats))';   // durees de sejour
traj = [etats;durees];

function traj = saut_harmonisee(L,x0,n)
//     Simule par la chaine harmonisee, n sauts du 
//     processus de generateur L a partir de x0. 
//     Retourne une matrice a n+1 colonnes et 2 lignes. 
//     La premiere ligne contient les etats visites, la
//     seconde contient les durees de sejour.

tm = max(-diag(L));              // horloge interne 
P = (L+tm*eye(L))./(tm*ones(L)); // matrice de transition
etats = grand(n,"markov",P,x0);  // etats visites
etats = [x0,etats];
durees = grand(etats,"exp",1);   // exponentielles
durees = durees/tm;              // durees de sejour
traj = [etats;durees];


