How does the high school approach to classifying linear systems generalize to more variables, and why does it work? If we have two linear equations in two variables, there may be $0$, $1$ or infinitely-many solutions. The standard approach to working out which case we're dealing with is probably to use Gaussian elimination. However, in high schools around here, they don't teach Gaussian elimination, and instead expect students to use (something like) the following result:

Theorem.
Consider $a,b,c,d,p,q \in \mathbb{R}$.
Let $\Gamma(x,y)$ denote the following condition. $$ax+by = p$$ $$cx+dy = q$$ 
Define $\tilde{\Gamma} := \{(x,y) \in \mathbb{R}^2 : \Gamma(x,y)\}.$
Then:
  
  
*
  
*$\tilde{\Gamma}$ has exactly one element if and only if  $ad\neq bc$.
  
*$\tilde{\Gamma}$ has infinitely many elements if and only if $ad = bc$ and $bq = dp.$
  
*$\tilde{\Gamma}$ has no elements if and only if $ad = bc$ and $bq \neq dp.$
(*I haven't actually proved this, and some nonzeroness assumptions might be necessary to make it fully technically correct.)

What seems to be happening, in particular, is that minors play a role in the number of solutions. By a $k$-minor of a matrix $M$, I mean a determinant of a $k \times k$ submatrix of $M$. What seems to be going on is that we have to look at the $2$-minors of the "augmented matrix"
$$
  \left[\begin{array}{rr|r}
    a & b & p \\
    c & d & q \\
  \end{array}\right]
$$
or, perhaps more precisely, the $2$-minors of 
$$
  \left[\begin{array}{rrr}
    a & b & -p \\
    c & d & -q \\
  \end{array}\right]
$$
which, due to multilinearity of the determinant in the columns of the matrices to which it can be applied, has the same $2$-minors up to a change in sign.
However, I don't really get what's going on here. I get that we can move to projective space by homogenizing the linear equations above to obtain $$ax+by -pz = 0$$ $$cx+dy-qz=0$$
This seems to have something do with it. I also get that there's a way to determine the rank of a matrix using its minors, and thus the nullity using the rank-nullity theorem. But this isn't enough to really understand what's going on.

Question. How does this classification generalize to more variables and/or more equations, and does projective space have something to do with it? Why does it work, and what's really going on here?

I put in some partially-related tags to get input from a wider audience, but feel free to remove them if you think they're a bad fit.
 A: So you're trying to classify the number of solutions of $MX = P$ where $P$ is a prescribed vector, $M$ a prescribed matrix, and $X$ is the unknown. 
Obviously the first condition you get is that there is a unique solution if and only if $\det(M) \neq 0$ (or "is invertible" if you're working over more general rings - though what I'll say later doesn't work on arbitrary rings as stated). Indeed, if $\det(M)=0$, then either $P$ is in the image in which case there is a whole affine space of solutions, or $P$ isn't, in which case there aren't solutions. In this situation, an explicit solution is given by the explicit formula for the inverse in terms of the determinant and the minors of $M$.
Now let's move on to $\det (M)= 0$. As I said, here, there are two cases, depending on whether $P$ is in the image of $M$ or not; and your goal would be to classify this with some minors.
But now note the following thing : the columns of $M$, say $C_1,...,C_k$ are precisely the image of the basis vectors, $e_1,...,e_k$, so the image of $M$ is the span of the columns. Therefore $P$ is in the image if and only if there is a linearly independent subfamily $C_{i_1},...,C_{i_m}$ such that $C_{i_1},...,C_{i_m},P$ is not linearly independent. 
In other words, you want there to be an extracted determinant of $M$ of size, say $m$, that isn't $0$, but such that the extracted determinant of size $m+1$ with the same columns, and any new set of rows and  adding $P$ (the corresponding rows of $P$) is $0$ (if you're not adding the condition on the determinant of size $m$, this doesn't tell you anything, because the columns might just be dependent). 
This gives, even in dimension $2$, a rather more complicated description of the solutions (first of all it won't only depend on $b,d$, there will also be a case analysis on $a,c$, and you will indeed have some questions about who is zero and who isn't). 
So the global plan of description with such a system of equations is : 
1- Compute the determinant, see if it's $0$ or not. If it's nonzero, you're done and you can extract solutions explicitly with the minors and the determinant. If it's $0$, move on to step 2.
2- For $m$ going from $1$ to $k-1$, compute extracted determinants of size $m$ (so choosing specific rows of a stystem of $m$ columns) and for the ones that aren't $0$, choose the same columns and add $P$  and see if it's $0$ with any choice of rows . 
3- If none of the linearly independent families of columns become dependent after adding $P$, then $P$ is not in the image and there are $0$ solutions.
4- If one of them does, say with the columns $C_{i_1},...,C_{i_m}$, then $P$ is in the image, more specifically it's in the image of those columns. With a bit of luck (or some more computation) you can find a specific antecedent $Y$, and then the set of antecedents if $Y+\ker M$.
So it has to do with the minors of $[ M \mid P ]$ indeed, because you want to see what families of columns of $M$ become dependent after adding $P$, which is measured by $(k-1)$- and $k$-minors, for $k$ varying a bit. 
I don't know enough about the projective plane to say if it jas anything to do with it.
A: I think the generalization of the method described in the question for the 2 dimensional case is just the following:


*

*Find the largest set of free vectors extracted from the columns of $f$, by computing minors. This gives the dimension of the solution set if it's not empty ($\DeclareMathOperator{\rk}{rk}\dim \ker(f) = n - \rk(f)$ with $n$ the dimension).

*By computing minors, test if adding $y$ to this set of vectors gives a free family of vectors. There is a solution if the answer is "no".


This can be simplified in some cases and becomes in the 2d case:


*

*We first remove the case where $f=0$. Thus $\rk(f) ∈ \{1,2\}$.

*If $\rk(f)=2$, then there is a unique solution (no need for step 2).

*If $\rk(f)=1$, then either $(a,c)≠(0,0)$ or $(b,d)≠0$. In the second case, there is a solution if and only if $bq≠dp$ (there is only one minor to compute in step 2).



Here is an explanation of the link with projective spaces.
$\newcommand{\bm}[1]{\begin{bmatrix}#1\end{bmatrix}}$
Let $f : V→W$ be a linear map and let $y ∈ W$. We want to solve $f(x)-y=0$. A vector $x∈V$ is a solution if and only if $\bm{f&-y}\bm{x\\1}=0$.
Note. The free vector space over the affine space $V$ is $V⊕ℝ$, in the sense of category theory. Under this viewpoint (adjunction), the affine map $x ∈ V ↦ f(x)-y ∈ W$ corresponds to the linear map $(x,p) ∈ V⊕ℝ ↦ f(x)-py ∈ W$ denoted above by $\bm{f&-y}$.
The question is to know the dimension of the intersection of $\ker(\bm{f&-y})$ and the set of vectors of the form $\scriptsize\bm{x\\1}$. The dimension of $\ker(\bm{f&-y})$ is $(n+1)-\DeclareMathOperator{\rk}{rk}\rk(\bm{f&-y})$, and $\rk(\bm{f&-y})$ is the size of the largest non-zero minor of $\bm{f&-y}$.
From a projective viewpoint, the dimension of the space of solutions is thus always $n-\rk(\bm{f&-y})$ (we remove $1$). The solutions "at infinity" are the kernel of $f$.
From an affine viewpoint, there are two possibilities: either there is no solution, or there is a solution, in which case the dimension of the space of solutions is the same as in the projective case. The second case arises if and only if $\rk(f) = \rk(\bm{f&-y})$.
