A = [3 0 3; 0 -1 3 ; 1 3 0];
b = [ 12 7 7 ]';

[m,n] = size(A);
if m ~= n 
    error('Matrix A must be square');
end

% Augmented coefficient matrix for equivalence transformations
Ab = [A b];

% Forward elimination
for i = 1: n-1                          % all rows treated
    for j = i+1:n                       % all elements below diagonal set 0
        l = Ab(j,i) / Ab(i,i);          % No pivoting
        Ab(j, :) = Ab(j,:) - l*Ab(i,:);
    end
end

% Backward substitution
x(n) = Ab(n,n+1) / Ab(n,n);

for i = n-1: -1 : 1
    x(i) = (Ab(i, n+1)- Ab(i, i+1:n)*x(i+1:n)) / Ab(i,i); 
end

x