A=[1, 1+0.5e-15, 3; 2,2,20; 3,6,4];
b=[5+0.5e-15;24;13];

% we solve by LU without pivoting
% [L,U]=lufact(A);
% y=L\b;
% x=U\y
% 
% xex=ones(3,1);
% 
% err=norm(xex-x)/norm(xex);

% we solve by LU with pivoting
[L,U,P]=lu(A);
y=L\(P*b);
x=U\y

xex=ones(3,1);

err=norm(xex-x)/norm(xex)