Finding a basis for Kernel and Image of a linear transformation using Gaussian elimination. This is a very fast method for computing the kernel/nullspace and image/column space  of a matrix. I learned this from my linear algebra teacher but I haven't seen it mentioned online apart from this reference in wikipedia.
The method is as follows: If $A$ is a m × n matrix, we construct $\displaystyle\left[\frac{I}{A}\right]$, where $I$ is the n × n identity matrix.
We then do elementary column operations until our $A$ is in column echelon form and we get the augmented matrix $\displaystyle\left[\frac{C}{B}\right]$. A basis of the kernel of $A$ consists in the columns of $C$ such that the corresponding column of $B$ is a zero column. A basis for the image of $A$ consist of all the non-zero columns of $B$.
An example:
$A =
               \begin{pmatrix}
               4 & 1  & 3 \\
               2 & -1  & 3 \\
                2 & 1 &  1 \\
                 1 & 1 & 0
                  \end{pmatrix}$
$\displaystyle\left[\frac{I}{A}\right] =  \begin{pmatrix}
               1 & 0 & 0 \\
               0 & 1 & 0 \\
               0 & 0 & 1 \\
               \hline 4 & 1  & 3 \\
               2 & -1  & 3 \\
                2 & 1 &  1 \\
                 1 & 1 & 0
                  \end{pmatrix}$ $\to$  $\begin{pmatrix}
               1 & 0 & -1 \\
               0 & 1 & 0 \\
               0 & 0 & 1 \\
               \hline 4 & 1  & -1 \\
               2 & -1  & 1 \\
                2 & 1 &  -1 \\
                 1 & 1 & -1
                  \end{pmatrix}$$\to$ $\begin{pmatrix}
               1 & 0 & -1 \\
               0 & 1 & 1 \\
               0 & 0 & 1 \\
               \hline 4 & 1  & 0 \\
               2 & -1  & 0 \\
                2 & 1 &  0 \\
                 1 & 1 & 0
                  \end{pmatrix}$
So a basis for $ Ker(A) = \begin{pmatrix} -1 \\ 1 \\ 1 \end{pmatrix}$ 
And a basis for $ Im(A) = \begin{pmatrix} 4 \\ 2 \\ 2 \\ 1 \end{pmatrix},\begin{pmatrix} 1 \\ -1 \\ 1 \\ 1 \end{pmatrix}$.
My problem is that I don't understand fully why this method works. From the wikipedia article I learned this derives from Gaussian Elimination, which sort of makes sense to me as it is similar to the method for finding the inverse of a matrix using Gauss-Jordan Elimination. So my question is: Why does this work?

The image of $A$ is precisely the column space of $A$, that is, the vector space spanned by the columns of $A$. Performing column operations does not change the column space of $A$ (see below), so the nonzero columns in the echelon form are a basis for the column space of $A$ (and, hence, for the image of $A$). 
A vector $(b_1,\dots,b_n)$ is in the nullspace of $A$ precisely when $b_1C_1+\cdots+b_nC_n=0$, where $C_1,\dots,C_n$ are the columns of $A$. The column operations you are doing are forming linear combinations of the columns of $A$, and the top part of your augmented matrix is keeping track of the coefficients of those linear combinations. So when you get a zero column, the entries above it are the coefficients of a vanishing linear combination of the columns of $A$; they form an element of the nullspace. 
"Performing column operations does not change the column space of $A$." There are three kinds of elementary column operation. You just have to check that for each of these three types, performing it does not change the column space. 
