%iterative method for
%mean first passage time tensor mu of an ergodic higher order chain
%P is m-th order, n-dimensional transition tensor
%mu0 is initial approximation to mu
%tol is tolerance
%mu0 and tol are optional, default mu0=ones(s), default tol=1e-6
%specify both mu0 and tol if any of them is different from its default
function mu=mfpti(P,mu0,tol)
s=size(P);
n=s(1);
m=length(s);
N=n^m;
if nargin<2
    mu0=ones(s);
    tol=1e-6;
end
kmax=1000;
for k=1:kmax
    mu=ones(s)+bprod(mu0-diagt(mu0),P);
    if norm(reshape(mu,1,N)-reshape(mu0,1,N),inf)<tol
        fprintf('mfpt algorithm converges in %3d iterations\n',k)
        break
    else
        mu0=mu;
    end
end
