function instab(epsilon)
n = 3
xvera = ones(n,1);
A = [1 1 -3; 2 2-epsilon 4; 1 9 4];
b = A* xvera;
[L,U] = lugauss(A);
disp('A - LU')
ERR = A - L*U
y = L\b;
x = U\y;
err_rel = norm(x-xvera)/norm(xvera);
res_rel = norm(b-A*x)/norm(b);

disp('senza pivot, x = ')
fprintf('\t %16.16f \n',x)
fprintf('\n\t err_rel= %12.12e res_rel= %12.12e \n',err_rel, res_rel)

[L,U,P] = lu(A);
disp('PA - LU')
ERR = P*A - L*U
y = L\(P*b);
x = U\y;
err_rel = norm(x-xvera)/norm(xvera);
res_rel = norm(b-A*x)/norm(b);


disp('con pivot, x = ')
fprintf('\t %16.16f \n',x)
fprintf('\n\t err_rel= %12.12e res_rel= %12.12e \n',err_rel, res_rel)
