function [x, Ainv] = Gauss_Jordan(A, b)
N = length(b);

B  = [A,  diag(ones(N,1)), b];

for n = 1:N;
    B(n,:) = B(n, :)./B(n,n);
    for m = 1: N
        if n ~= m
            B(m,:) = B(m,:) - B(n,:)*B(m,n);
        end;
    end;
    disp(B);
end;
x = B(:,2*N+1);
Ainv = B(:, (N+1):(N+N));



