Help with converting indexed computation to matrix multiply

Hi there, I'm pretty stuck on converting an indexed computation to a matrix multiply.

I currently have:

Pkjn = zeros(N,1);
A = zeros(N,N);
s = zeros(3,N);
W = zeros(3,N);
n = 0:N-1;
for k=1:N
    for j=1:N               
        Pkjn = cos((k-j)*2*pi*n/N).';  

          Mk = 3 * s(:,k); %3x1
          Mj = 3 * s(:,j); %3x1

          A(k,j) = (repmat(Mj.*Mk,1,N).*W  * Pkjn ).' * ones(3,1); 
      end
  end

and I'd like to convert the computation of A to a Matrix multiply operation. I've got as far as:

h = dftmtx(N);
a = ctranspose(h)*h;

for cc=1:3
    b(:,:,cc)=9*s(cc,:)'*(s(cc,:)).*a;
end
A = sum(b,3)

but I can't seem to get any further. Any helpful suggestions gratefully received!

Thanks, Ed.