A=[4 0 -8; 1 2 3; -2 4 -1]

% switch rows 1 and 2
P=[0 1 0; 1 0 0; 0 0 1]
P*A

A([1 2],:)=A([2 1],:)

%-------------------------------------------------------
% 4-component system where pivoting is necessary
A=[1 2 3 -2; 2 4 1 0; 3 3 2 5; -1 6 2 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,:)
 
L(4,1)=U(4,1)/U(1,1)
U(4,:)=U(4,:)-L(4,1)*U(1,:)

%-------------------------------------------------------
% first non-zero pivot: switch rows 2 and 3
U([3 2],:)=U([2 3],:)  % switch rows 2 and 3 in U
L([3 2],1)=L([2 3],1)  % switch the corresponding saved multipliers in L

%P=[1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]  % permutation matrix
P=eye(size(A))
P([3 2],:)=P([2 3],:)

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

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

% check that L*U=P*A
L*U-P*A

%-------------------------------------------------------
% maximal entry pivot: switch rows 2 and 4
U([4 2], :)=U([2 4],:)  % switch rows 2 and 4 in U
L([4 2], 1)=L([2 4],1)  % switch the corresponding saved multipliers in L

P=eye(size(A)); P([4 2],:)=P([2 4],:) % permutation matrix

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

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

% check that L*U=P*A
L*U-P*A
