Why solving a system of linear equation produces the intersection of the equation 1) $x+y=1$
2) $-x+y=1$
Geometrically we can visualize the two lines will intersect at $x=0, y=1$. Consider this algebraic solution using Gaussian Elimination,
.
But why do they be reduced to the intersection point?
Is there any algebraic mathematical proof for solution of linear equations?
 A: Let the lines be given by
$$
y=1-x;\quad y=1+x.\tag{*}
$$
If a pair $(x_0,y_0)$ satisfies
$$
x_0+y_0=1,\quad -x_0+y_0=1
$$ 
then $(x_0,y_0)$ simultaneously satisfies both the equations in (*), so the pair represents a point on the $(x,y)$-plane that is an intersection between the 2 lines.
A: So you start with $$x+y=1$$$$-x+y=1$$
Keep the first equation the same, and add the first equation to the second equation to get a new second equation: $$x+y=1$$$$ 2y=2$$
These hold because the original equations hold, and the step is reversible, so they are equivalent to the original pair and have the same solutions. Now divide the second equation by $2$ (also reversible) to obtain $$x+y=1$$$$ y=1$$ Now subtract the second equation from the first $$x=0$$$$ y=1$$
By taking reversible steps we have shown that the final solution is equivalent to the original equations.
This Gaussian elimination reflects the steps you would take to solve the equation. It is simply another way of writing the same thing. The advantage is that it does this in a controlled and systematic manner which applies to all cases - we can write an algorithm for it. And this makes it possible to analyse the process mathematically, and thus to abbreviate it as a matrix solution, and to analyse cases where no solution is possible.
A: For a proof of how the Gaussian Elimination works take a look at Linear Algebra - Hoffman, Kunze. Reading Chapter 1 and Section 2.5 would do. It is rather lengthy but the gist is this. 
We are allowed three elementary row operations. Non-zero scalar multiplication, interchanging rows and adding a scalar multiple of one row to the other. The resulting system after one of these operations are applied is a "linear combination of rows from the first system". So any solution to the first is also a solution to the new system (see the reference for more details). Now, the row operations are defined in such a way that they all have an inverse row operation - respectively, multiplying by the inverse of the aforementioned non-zero scalar, interchanging of the same two rows and adding the negative time the aforementioned scalar multiple of the row to the second row. Now these inverse operations restore the old system and they too are "linear" operations on the rows. It can be visualised then how all solutions to the new system are solutions to the original one. Hence, the solutions to two systems are the same if one is obtained by applying one (or a finite number of ) row operation/s to the other. We reduce a matrix to an Echelon matrix so we can read the solution off easily. As is the case with reducing to an identity matrix. 

Now you have above a $2 \times 2$ system. So the equations can be represented on the plane. Note that a solution to the system is a pair $(x, y)$ such that both equations are satisfied. The line on the plane given by the equation $ x + y = 1 $ represents all points on $\Bbb R \times \Bbb R$ which satisfy the equation. Same for the second equation. Now a solution to the system must satisfy both equations. So geometrically it should be a point on both lines. If the lines are not overlapping or are parallel there can be only one such point which is the intersection point. 
How the identity matrix gives the solution is explained roughly at the top. But you will have to read more to get a rigorous explanation. The book recommended here is excellent for the purpose. 
