You won't need matrices of different size.
I suggest to use 4x3 matrices which is a 3x3 matrix with an additional translation-vector.
Most APIs use 4x4 matrices, but I'm pretty sure you won't need those for a while.
Matrix:
xx xy xz xp
yx yy yz yp
zx zy zz zp
[rotation] [position]
Try to understand that the 3x3 part consists of three vectors.
These vectors are the axes of your coordinate system:
(xx,yx,zx) points to the right
(xy,yy,zy) points up
(xz,yz,zz) points into depth
Identity matrix (transforms x,y,z to the same x,y,z):
1 0 0 0
0 1 0 0
0 0 1 0
Notice that the bases are the cartesian axes (transformation to the vector-space you're already in).
A vector describes a position or direction relative to a coordinate system.
When transforming a vector (x,y,z) to another coordinate-system, the components of the vector are then relative to new coordinate axes:
x'= xx*x + xy*y + xz*z + xp
y'= yx*x + yy*y + yz*z + yp
z'= zx*x + zy*y + zz*z + zp
Now we want to transform using two matrices "m1" and "m2":
first transform (x,y,z) to (x1,y1,z1) using m1 (using equation above)
x1= m1.xx*x + m1.xy*y + m1.xz*z + m1.xp
y1= m1.yx*x + m1.yy*y + m1.yz*z + m1.yp
z1= m1.zx*x + m1.zy*y + m1.zz*z + m1.zp
then transform (x1,y1,z1) to (x2,y2,z2) using m2
x2= m2.xx*x1 + m2.xy*y + m2.xz*z + m2.xp
y2= m2.yx*x1 + m2.yy*y + m2.yz*z + m2.yp
z2= m2.zx*x1 + m2.zy*y + m2.zz*z + m2.zp
This is twice the work, so we're looking for new matrix "n" which does both transforms at once.
We simply plug the first equation into the second and reorder back to matrix-layout (just shown for one component):
x2= m2.xx * (m1.xx*x + m1.xy*y + m1.xz*z + m1.xp) (equation for x1)
+ m2.xy * (m1.yx*x + m1.yy*y + m1.yz*z + m1.yp) (equation for y1)
+ m2.xz * (m1.zx*x + m1.zy*y + m1.zz*z + m1.zp) (equation for z1)
+ m2.xp
=> remove braces
x2= m2.xx*m1.xx*x + m2.xy*m1.yx*x + m2.xz*m1.zx*x
+ m2.xx*m1.xy*y + m2.xy*m1.yy*y + m2.xz*m1.zy*y
+ m2.xx*m1.xz*z + m2.xy*m1.yz*z + m2.xz*m1.zz*z
+ m2.xx*m1.xp + m2.xy*m1.yp + m2.xz*m1.zp + m2.xp
=> refactor to x,y,z:
x2= (m2.xx*m1.xx+m2.xy*m1.yx+m2.xz*m1.zx) * x
+ (m2.xx*m1.xy+m2.xy*m1.yy+m2.xz*m1.zy) * y
+ (m2.xx*m1.xz+m2.xy*m1.yz+m2.xz*m1.zz) * z
+ (m2.xx*m1.xp + m2.xy*m1.yp + m2.xz*m1.zp + m2.xp)
So here we've got the first row of a matrix which transform "x" directly to "x2".
The complete new matrix "n" looks like this:
n.xx= m2.xx*m1.xx + m2.xy*m1.yx + m2.xz*m1.zx
n.xy= m2.xx*m1.xy + m2.xy*m1.yy + m2.xz*m1.zy
n.xz= m2.xx*m1.xz + m2.xy*m1.yz + m2.xz*m1.zz
n.xp= m2.xx*m1.xp + m2.xy*m1.yp + m2.xz*m1.zp + m2.xp
n.yx= m2.yx*m1.xx + m2.yy*m1.yx + m2.yz*m1.zx
n.yy= m2.yx*m1.xy + m2.yy*m1.yy + m2.yz*m1.zy
n.yz= m2.yx*m1.xz + m2.yy*m1.yz + m2.yz*m1.zz
n.yp= m2.yx*m1.xp + m2.yy*m1.yp + m2.yz*m1.zp + m2.yp
n.zx= m2.zx*m1.xx + m2.zy*m1.yx + m2.zz*m1.zx
n.zy= m2.zx*m1.xy + m2.zz*m1.zy + m2.zy*m1.yy
n.zz= m2.zx*m1.xz + m2.zy*m1.yz + m2.zz*m1.zz
n.zp= m2.zx*m1.xp + m2.zy*m1.yp + m2.zz*m1.zp + m2.zp