% function [x,LU]=Gausseli(A,b) - Gaussian Elimination without par pivoting
%
function [x,LU]=Gausseli(A,b)
% A=[ 2,  1, -2; -4, -1, 5; 2, 2, 2];
% b=[-2;  3;  0];
% x=[ 1; -2;  1]; Solution for A*x=b
[m n]=size(A);
if (m~=n)
    error('matrix A must be square');
end
nb=n+1;
Aug=[A, b];
% Forward Elimination
for i=1:n-1
    for k=i+1:n
        r=Aug(k,i)/Aug(i,i);
        Aug(k,i+1:nb)=Aug(k,i+1:nb)-r*Aug(i,i+1:nb);
        Aug(k,i)=r;
    end
end
% Back Substitution
x=zeros(n,1);
x(n)=Aug(n,nb)/Aug(n,n);
for j=n-1:-1:1
    x(j)=(Aug(j,nb)-Aug(j,j+1:n)*x(j+1:n))/Aug(j,j);
end 
LU=Aug(1:m,1:n);
% Output x=[1, -2, 1]';
% LU=L*U in matrix Aug(1:m, 1:n)
%  2  1 -2             1  0  0        2  1 -2
% -2  1  1, where L = -2  1  0 ;  U = 0  1  1
%  1  1  3             1  1  1        0  0  3