%block diagonal direct method for
%mean first passage time tensor mu of an ergodic higher order chain
%P is m-th order, n-dimensional transition tensor
function mu=mfptd(P)
s=size(P);
mu=zeros(s);
n=s(1);
m=length(s);
N1=n^(m-1);
b=ones(N1,1);
if m>2
    s=s(2:end);
else
    s=[s(2:end) 1];
end
ind=cell(1,m-1);
for k=1:n
    C=eye(N1);
    for i=1:N1
        [ind{:}]=ind2sub(s,i);
        indi=fliplr([ind{:}]);
        for j=1:n
            if j~=k
                indP=num2cell([j indi]);
                indc=num2cell(fliplr([j indi(1:end-1)]));
                indc=sub2ind(s,indc{:});
                C(i,indc)=C(i,indc)-P(indP{:});
            end
        end
    end
    muk=C\b;
    for i=1:N1
        [ind{:}]=ind2sub(s,i);
        indi=fliplr([ind{:}]);
        indmu=num2cell([k indi]);
        mu(indmu{:})=muk(i);
    end
end
        
        
