Here's a recipe for matrix multiplication that you can't forget: choose bases b_i/c_j for your domain/codomain. Then all a matrix is is listing the outputs of a function for your basis: if you have a linear function f, then the ith column of its matrix A is just f(b_i). If you have another function g from f's codomain, then same thing, its matrix B is just the list of outputs g(c_j). Then the ith column of BA is just g(f(b_i)). If you write these things down on paper and expand out what I wrote, you'll see the usual row and column thing pop out. The point is that f(b_i) is a weighted sum of the c_i (since c_i is a basis for the target of f), but you can pull the weighted sums through the definition of g because it's linear. A basis gives you a minimal description/set of points where you need to define a function, and the definition for all other points follows from linearity.

The point of the eigen-stuff is that along some directions, linear functions are just scalar multiplication: f(v) = av. If the action in a direction is multiplication by a, then it can't also be multiplication by b. So unequal eigenvalues must mean different directions/linearly independent subspaces. So e.g. if you can find n different eigenvalues/eigenvectors, you've found a simple basis where each direction is just multiplication. You also know that it's invertible if the eigenvalues are nonzero since all you did was multiply by a_i along each direction, so you can invert it by multiplying by 1/a_i on each direction.

Taught properly it's all very straightforward, though determinants require some more buildup with a detour through things like quotienting and wedge products if you really want it to be straightforward IMO. You start by saying you want to look at oriented areas/volumes, and look at the properties you need. Then quotienting gives you a standard tool to say "I want exactly the thing that has those properties" (wedge products). Then the action on wedges gives you what your map does to volumes, with the determinant as the action on the full space. You basically define it to be what you want, and then you can calculate it by linearity/functoriality just like you expand out the definition of a linear map from a basis.