function X = gauss_jordan(A,B) 

% function X = gauss_jordan(A,B) 
% solution of AX = B

X = [A, B];
for k = 1:size(A,1);
   [xmax, imax] = max(abs(X(k:end,1)));
   imax = imax(1)+k-1;
   X([k, imax],:) = X([imax, k],:);
   v = X(:,1); v(k) = 0;
   X(k,:) = X(k,:)/X(k,1);
   X = X(:,2:end) - v*X(k,2:end);
end;

