
function [exact_sol,num_sol] = gauss_back(A,b);

exact_sol = inv(A)*b;

% gauss elimination and backward substitution

A = [A,b];

[n,m] = size(A);
TOL = 1.0e-8;

if (n ~= m-1)
    display('the size of A is not correct\n');
end

for i = 1:n-1;
    
% find the maximum value for i-th step   
    p = i;
    max_val = abs(A(i,i));
    for k = i+1:n;
        if (abs(A(k,i)) > max_val)
            p = k;
            max_val = abs(A(k,i)); 
        end
    end
    
    %[max_val,p] = max(abs(A(i:n,i)));
    %p = p+i-1;

% switch two rows
    
    if (p ~= i)
        temp = A(i,:);
        A(i,:) = A(p,:);
        A(p,:) = temp;
    end
    
    if (abs(A(i,i)) < TOL)
        display('The solution is not unique\n');
        return;
    end
    
    for j = (i+1):n;
        m = A(j,i)/A(i,i);
       A(j,:) = A(j,:)-m*A(i,:);
    end
end

x(n) = A(n,n+1)/A(n,n);

for i = (n-1):-1:1;
    SUM = 0;
    for j = (i+1):n;
        SUM = SUM+A(i,j)*x(j);
    end
    x(i) = (A(i,n+1)-SUM)/A(i,i);
end
num_sol = x';
