%Linear equation system 'Ax=r' by Gauss elimination method.
%Written by: "Sobhan Rostami"
%MSc student of structure engineering of Azad university of kerman.
clc
clear all
%=================================================================
disp('Solution of N-equation "[A][X]=[r]"')
A   = [2 2 -1 1; -1 1 2 3; 3 -1 4 -1; 1 4 -2 2]
b   = [7 3 31 2]
r   = b'
n   = length(A)
%-----------------------------------------------------------------
%create upper triangular matrix
s=0;
for j=1:n
    if A(j,j)==0
        k=j;
        for k=k+1:n
            if A(k,j)==0
                continue
            end
            break
        end
        B=A(j,:); C=r(j)
        A(j,:)=A(k,:); r(j)=r(k);
        A(k,:)=B; r(k)=C;
    end
    for i=1+s:n-1
        L=A(i+1,j)/A(j,j);
        A(i+1,:)=A(i+1,:)-L*A(j,:);
        r(i+1)=r(i+1)-L*r(j);
    end
    s=s+1;
end
%-----------------------------------------------------------------
%Solution of equations
x(n)=r(n)/A(n,n);
for i=n-1:-1:1
    sum=0;
    for j=i+1:n
        sum=sum+A(i,j)*x(j);
    end
    x(i)=(1/A(i,i))*(r(i)-sum);
end
%------------------------------
%Checking with matlab functions
p=inv(A)*r;
disp(A)
%------------------------------
%Output
disp('@----------------------------------------------------------@')
disp('Output [B][x]=[b]')
disp('Upper riangular Matrix [B] =');disp(A)
disp('Matrix [b] =');disp(r)
disp('solution of linear equations :');disp(x')
disp('solve with matlab functions(for checking):');disp(p)