function AB = tensprod(A,B)
[mA,nA] = size(A);
[mB,nB] = size(B);
AB = [];
for x=1:nA
  for y = 1:nB
    abind = coord2ind([x y],[nA nB]);
    AB(:,abind) = tensprod1(A(:,x),B(:,y));
  end
end
return


function pq = tensprod1(p,q)
mp = length(p);
mq = length(q);
pq = zeros(mp*mq,1);
for x = 1:mp
  for y = 1:mq
    xyind = coord2ind([x y],[mp mq]);
    pq(xyind) = p(x) * q(y);
  end
end
return