Creating a matrix from a dyadic product Ax, Ay, Az: [N-by-N]
B=AA (dyadic product)  


means  

B(i,j)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)] 


B(I,j): a 3x3 matrix.

One way to construct B is:

N=2;
Ax=rand(N); Ay=rand(N); Az=rand(N);
t=1;
F=zeros(3,3,N^2);
for  i=1:N
    for j=1:N
        F(:,:,t)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)];
        t=t+1;      %# t is just a counter
    end
end
%# then we can write
B = mat2cell(F,3,3,ones(N^2,1));
B = reshape(B,N,N)'; 
B = cell2mat(B);


Is there a faster way than this, especially when N is large?
There is a problem with vector multiplication in the second loop. You should transpose the second vector before doing multiplication.

F(:,:,t)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)]';


One of the ways to speed things up is to apply vectorization instead of two for loops, such as with Dyadics.