%direct method for computing the fundamental tensor sigma
%of an (m-1)-th order absorbing chain
%with absorbing states t+1, t+2, ..., n,
%whose transition tensor is P=[Q O; R I],
%Q is an m-th order, t x t x n x ... x n tensor,
%R is an m-th order, a x t x n x ... x n tensor,
%O is an m-th order, t x a x n x ... x n zero tensor,
%and I is an m-th order, a x a x n x ... x n identity tensor
%m>=2, n>=2, 1<=t<n, and a=n-t
function sigma=fund(Q)
s=size(Q);
sigma=zeros(s);
n1=s(1);
m=length(s);
if m>2
    s=fliplr(s(2:end));
else
    s=[s(2:end) 1];
end
N1=prod(s);
ind=cell(1,m-1);
for k=1:n1
    C=eye(N1);
    b=zeros(N1,1);
    for i=1:N1
        [ind{:}]=ind2sub(s,i);
        indi=fliplr([ind{:}]);
        if indi(1)==k
            b(i)=1;
        end
        for j=1:n1
            indQ=num2cell([j indi]);
            indc=num2cell(fliplr([j indi(1:end-1)]));
            indc=sub2ind(s,indc{:});
            C(i,indc)=C(i,indc)-Q(indQ{:});
        end
    end
    sigmak=C\b;
    for i=1:N1
        [ind{:}]=ind2sub(s,i);
        indi=fliplr([ind{:}]);
        indsigma=num2cell([k indi]);
        sigma(indsigma{:})=sigmak(i);
    end
end