%box product of m-th order, n-dimensional tensors A and B
function C=bprod(A,B)
s=size(A);
C=zeros(s);
n=s(1);
m=length(s);
N=n^m;
ind=cell(1,m);
for i=1:N
    [ind{:}]=ind2sub(s,i);
    indi=[ind{:}];
    indC=num2cell(indi);
    for j=1:n
        indA=num2cell([indi(1) j indi(2:end-1)]);
        indB=num2cell([j indi(2:end)]);
        C(indC{:})=C(indC{:})+A(indA{:})*B(indB{:});
    end
end


