I solved it myself
function result = mytensorprod(A,B,dims)
% performs the tensor product for matrices A and B
% A and B are contracted along the dimensions specified in dims
% First A and B are permuted to perform the inner product
% Same functionality as the matlab function tensorprod
sizeA = size(A); sizeB = size(B);
dimAouter = 1:length(sizeA); dimAouter(dims)=[];
dimBouter = 1:length(sizeB); dimBouter(dims)=[];
Ap = permute(A,[dimAouter,dims]);
Apr = reshape(Ap,[prod(sizeA(dimAouter)),prod(sizeA(dims))]);
Bp = permute(B,[dims,dimBouter]);
Bpr = reshape(Bp,[prod(sizeB(dims)),prod(sizeB(dimBouter))]);
result = Apr*Bpr;
result = reshape(result,[sizeA(dimAouter),sizeB(dimBouter)]);
end
