What it is equation of plane in three dimensional space of that most fits the points
A (0,0,1)
B (1,0,1)
C (2,-1,0)
D (-1,1,0)
First, let's set up the system using matrices:
$$A = \begin{pmatrix} 0&0&1\\ 1&0&1\\ 2&-1&1\\ -1&1&1 \end{pmatrix}, X = \begin{pmatrix} x\\ y\\ c \end{pmatrix}, Y = \begin{pmatrix} 1\\ 1\\ 0\\ 0 \end{pmatrix}$$ Such that $$AX \approx Y$$
Solving using Least Squares: $$AX \approx Y$$ $$A^TAX=A^TY$$ $$X=(A^TA)^{-1}A^TY$$ $$X=\begin{pmatrix} 0\\ 0\\ 1/2 \end{pmatrix}$$
This tells us $0x+0y+1/2=z$, or $z=1/2$ is the "plane of best fit."
Note: please comment if I made a calculation error.
If you're looking for the best-fitting plane, using a least-squares estimate is one way to go.
Since the equation of a plane is $ax + by + cz = 1$, create a matrix $A$ where each row of $A$ is $[x_1\,y_1\,z_1]$. Since you have 4 data points, $A$ will have 4 rows.
Then, find the least squares solution to $Ax = 1$ where $x = [a\,b\,c\,]^T$.
So in your case,
$$A = \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 1 \\ 2 & -1 & 0 \\ -1 & 1 & 0 \\ \end{bmatrix}$$
First, three points are enough to define a plane, so you may ignore one of them for now. Plane equation is $<\vec{x},\vec{n}> = c$. $\vec{n}$ is a vector normal to the plane. Find it by picking any two vectors on the plane and taking their cross product. Make sure that the vectors aren't colinear though, or you'll get $\vec{0}$. Once you have $\vec{n}$, you can find $c$ by substituting the coordinates of any point for $\vec{x}$. Finally, you should check if this equation is satisfied by the point you chose at the beginning (if it is not, then the points aren't in one plane).