Your setup code must not be correct.
Consider first i=2.
In the second asignment chain(1,2,:) is assigned and has to do with the scalar "i" (there would no point looping of none of the entries involved "i"). The assignment also has to do with the vector "k", which is a 20 element vector, so the right hand side produces a vector of length 20, that is assigned to chain(1,2,:); this makes chain to be at least 1 x 2 x 20. The next assignment, to chain(2,1,:) is likewise of length 20, so chain must be at least 2 x 2 x 20.
Let us assume for the moment that in your 4th assignment you meant the left-hand side to be chain(2,2,:) instead of overwriting chain(1,2,:) that you assigned 2 lines above.
Now, go on to the next value of "i". The second assignment, to chain(1,2,:), has to do with "i" (changed) and the vector k (unchanged), producing a vector of length 20. But look where it is being assigned -- to chain(1,2,1:20), which was already assigned to for the first "i" !
Examining further, the same problem applies to all four assignments, so the effect of the code is as if it had only been executed for i=N (the last value).
As you said "Terms A,B,C and D of the matrix are of size 2x2x500", even though they individually are not 2 x 2 by anything, we can speculate that what you want is for "chain" to end up 2 x 2 x 500. "k" is length 20, and you have 50 distinct N values, 2 to 51, and 20 * 50 = 1000, exactly twice as long as desired. So maybe you meant that "chain" should come out 2 x 2 x 1000 ?
If that is what you want, then assign to chain(1,1,:,i-1) and chain(1,2,:,i1) and so on, producing a 2 x 2 x 20 x 50 array. You can then reshape() that to be 2 x 2 x 1000 and do the multiplications based upon that.
I will not advise on the multiplication at this time. Algebraic matrix multiplication is order-dependent, so it is crucial that you specify the order the matrices are to be multiplied in, and whether the multiplications are left-multiplications or right-multiplications. Or perhaps you want element-by-element multiplication instead of matrix multiplication? Seems less likely to me as it appears you are constructing rotation matrices.
