X=[1 -1 -1; 1 1 -1; 1 -1 1; 1 1 1];
disp(X);
M=X'*X;
D=M^(-1);
Y=[1;2;3;4];
disp(Y);
B=D*X'*Y;
disp(B);

N=7
t=linspace(-1,1,N);

for i=1:N
    for j=1:N
        z(i,j) = B(1) + B(2) * t(i) + B(3) * t(j);
    end
end

contour2d(t,t,z,N);
