function [A b x]=gausselk(A,b)
  n=length(b);
  augm=[A b]
  for k=1:n-1
    for i=k+1:n
      m=A(i,k)/A(k,k)
      for j=k:n
	A(i,j)=A(i,j)-m*A(k,j);
      endfor
      b(i)=b(i)-m*b(k);
    augm=[A b]
    endfor
  endfor

x=inv(A)*b;

endfunction