%box product of m-th order u x w x n x ... x n tensor A and w x v x n x ... x n tensor B
%where v must not exceed n
function C=bprode(A,B)
sA=size(A);
sB=size(B);
s=[sA(1) sB(2) sA(3:end)];
C=zeros(s);
w=sA(2);
m=length(s);
N=prod(s);
ind=cell(1,m);
for i=1:N
    [ind{:}]=ind2sub(s,i);
    indi=[ind{:}];
    indC=num2cell(indi);
    for j=1:w
        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


