/* inputs a discrete-time Markov transition matrix t and returns the 
   steady-state probabilities vector pi
   
   the basic idea is

   given pi * t = pi
   so pi * (t - I) = 0 (I is the identity matrix)
   also pi * 1v = 1  (1v is a vector of all 1s)
   then pi * tM = (0,0,...,0,1) 
      where tM is t - I but modified so that the last column is
      all 1s 
*/

Function("solvedmc",{t})  [
   nrows:=Length(t); 
   ForEach(i,1 .. nrows-1) [t[i][i]:=t[i][i]-1;]; 
   ForEach(i,1 .. nrows) [t[i][nrows]:=1;]; 
   v:=0*(1 .. nrows); 
   v[nrows]:=1; 
   SolveMatrix(Transpose(t),v); ];