function [x] = gauss_elim(A,b)
%  Solves : Ax=b using gaussian elmination
%  Returns: vector x
 
n  = max(size(A));
% Transform A to upper triangle matrix
for j=2:n,
    for i=j:n,
       m = A(i,j-1)/A(j-1,j-1);
       A(i,:) = A(i,:) - A(j-1,:)*m;
       A
       b(i) = b(i) - b(j-1)*m;
    end
end

% Back substitution

x = zeros(n,1);
x(n) = b(n)/A(n,n);

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