function x = gausselim(A, b)

% first, get a triangular system

n = size(A, 1);

if n ~= size(A, 2)
    error('fix your A matrix!')
end

for i=1:n
    
    aii = A(i,i);
    A(i,:) = A(i,:)/aii;
    b(i) = b(i)/aii;
    
    for j=(i+1):n
        
        aji = A(j,i);
        A(j,:) = A(j,:) - aji*A(i,:);
        b(j) = b(j) - aji*b(i);
        
    end
    
end

A
b