The easiest way to describe matrix multiplication is nested function composition.

    f(x) = 2x
    g(x) = x + 5

It should hopefully be obvious that "nesting" the two isn't commutative:

    g(f(x)) = (2x) + 5 = 2x + 5
    f(g(x)) = 2(x + 5) = 2x + 10

One miraculous fact here is that no matter how many functions we stack, we only ever have two terms: x and a constant term. Thus, we can represent this 'linear system' in terms of its coefficients, as long as we agree on an order: e.g. 2x+5 might become "2,5" in our system.

You can do "multiplication" on these packs to compose them together, and even though the rules feel obtuse in abstract, they follow the logic of function composition. A matrix represents a similar function transform, only it's in 3 dimensions, and in order to handle rotation, it needs to swap around x/y/z. So an identity matrix is really saying:

   f(x) = 1x + 0y + 0z + 0
   f(y) = 0x + 1y + 0z + 0
   f(z) = 0x + 0y + 1z + 0