L1=[2 0 0;1 -1 0; 3 4 5];
b=[1 2 3]';
x = forwsub(L1,b)

L2=[2 0 0;1 0 0; 3 4 5];
x = forwsub(L2,b)

U=[1 2 2; 0 -4 -6; 0 0 -3];
v=[3 -6 1]';
x=backsub(U,v)

A=[1 2 2; 4 4 2; 2 6 4];
b=[3 6 10]';
x=inv(A)*b

A=[1 2 2; 4 4 2 ; 2 6 4];
b=[3 6 10]';
[U,v]=gausselim(A,b)

M1= [1 0 0; -4 1 0; -2 0 1];
L1=inv(M1)

M2=[1 0 0; 0 1 0; 0 0.5 1];
L2=inv(M2)

L=L1*L2

A=[1 2 2; 4 4 2; 2 6 4];
[L,U]=lufactor(A)

delta=1/4;
A=[delta 1; 1 eps];
[L,U]=lufactor(A);
(L*U==A)

delta=0.5*eps;
A=[delta 1; 1 eps];
[P,L,U]=plufactor(A)
(P*A==L*U)
