% step-by-step LU factorization

%-------------------------------------------------------
% 3-component system, first class example
A=[4 0 -8; 1 2 3; -2 4 -1]

% initialize L and U
L=eye(size(A))
U=A

% Gauss elim column 1
L(2,1)=U(2,1)/U(1,1)
U(2,:)=U(2,:)-L(2,1)*U(1,:)

L(3,1)=U(3,1)/U(1,1)
U(3,:)=U(3,:)-L(3,1)*U(1,:)
 
% Gauss elim column 2
L(3,2)=U(3,2)/U(2,2)
U(3,:)=U(3,:)-L(3,2)*U(2,:)

% check L*U is equal to A
L*U

% solve A*x=b by forward and backward substitution
b=[12;2;7]
%b=[-9;8;23]

y=forsub(L,b)
x=backsub(U,y)

%-------------------------------------------------------
% 3-component system illustrating round-off error
format long
A = [10 -7 0; -3 2.09999 6; 5 -1 5]

U=A;
L=eye(size(A));

% Gauss elim first column
L(2,1)=U(2,1)/U(1,1)
U(2,:)=U(2,:)-L(2,1)*U(1,:)

L(3,1)=U(3,1)/U(1,1)
U(3,:)=U(3,:)-L(3,1)*U(1,:)
 
% Gauss elim second column

L(3,2)=U(3,2)/U(2,2)
U(3,:)=U(3,:)-L(3,2)*U(2,:)

b=[7; 3.90001; 6]

y=forsub(L,b)
x=backsub(U,y)
