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 

