A = chemreaction();
L = [1 0 0 0; 2 1 0 0; 2 1/3 1 0; 0 1/3 -1/5 1];
U = [1 0 -1 0; 0 3 -4 -1; 0 0 10/3 -5/3; 0 0 0 0];

norm(L*U - A)

x = [1; 2; 1; 2];
norm(A*x)
