n = 4;
A1 = diag( 3 * ones( n, 1 ), 0 ) + diag( - 2 * ones( n - 1, 1 ), 1 ) ...
       + diag( - 1 * ones( n - 1, 1 ), - 1 );
A1_det_sub = []; 
for p = 1 : n - 1
    A1_det_sub = [ A1_det_sub, det( A1( 1 : p, 1 : p ) ) ];    
end
A1_det_sub
% A1_det_sub =
%    3.0000    7.0000   15.0000

A2 = hilb( n );
format short e
A2_eig = eig(A2)'
% A2_eig =
%     9.6702e-05   6.7383e-03   1.6914e-01   1.5002e+00
format 

n = 9;  % Matrix A_1
A1 = diag( 3 * ones( n, 1 ), 0 ) + diag( - 2 * ones( n - 1, 1 ), 1 ) ...
       + diag( - 1 * ones( n - 1, 1 ), - 1 );
x1_ex = ones( n, 1 );       b1 = A1 * x1_ex;
[ L1, U1, P1 ] = lu( A1 );
[ y1 ] = forward_substitutions( L1, P1 * b1 );
[ x1 ] = backward_substitutions( U1, y1 );
e1_rel = norm( x1_ex - x1 ) / norm( x1_ex )
% e1_rel =
%    1.2820e-16
r1_rel = norm( b1 - A1 * x1 ) / norm( b1 )
% r1_rel =
%    3.5804e-16
A1_cond = cond(A1)
% A1_cond =
%    26.3742

A2 = hilb( n );   % Matrix A_2
x2_ex = x1_ex;   b2 = A2 * x2_ex;
[ L2, U2, P2 ] = lu( A2 );
[ y2 ] = forward_substitutions( L2, P2 * b2 );
[ x2 ] = backward_substitutions( U2, y2 );
e2_rel = norm( x2_ex - x2 ) / norm( x2_ex )
% e2_rel =
%      1.0906e-05
r2_rel = norm( b2 - A2 * x2 ) / norm( b2 )
% r2_rel =
%      1.3242e-16
A2_cond = cond(A2)
% A2_cond =
%      4.9315e+11

c1_v = [ ];   e1_rel_v = [ ];   r1_rel_v = [ ];
c2_v = [ ];   e2_rel_v = [ ];   r2_rel_v = [ ];
n_vect = 4 : 13;
for n = n_vect
    A1 = diag( 3 * ones( n, 1 ), 0 ) + diag( - 2 * ones( n - 1, 1 ), 1 ) ...
      + diag( - 1 * ones( n - 1, 1 ), - 1 );
    x1_ex = ones( n, 1 );     b1 = A1 * x1_ex; 
    [ L1, U1, P1 ] = lu( A1 );
    [ y1 ] = forward_substitutions( L1, P1 * b1 );
    [ x1 ] = backward_substitutions( U1, y1 );
    e1_rel_v = [ e1_rel_v, norm( x1_ex - x1 ) / norm( x1_ex ) ];
    r1_rel_v = [ r1_rel_v, norm( b1 - A1 * x1 ) / norm( b1 ) ];
    c1_v = [ c1_v, cond( A1 ) ];
    
    A2 = hilb( n );
    x2_ex = x1_ex;    b2 = A2 * x2_ex;
    [ L2, U2, P2 ] = lu( A2 );
    [ y2 ] = forward_substitutions( L2, P2 * b2 );
    [ x2 ] = backward_substitutions( U2, y2 );
    e2_rel_v = [ e2_rel_v, norm( x2_ex - x2 ) / norm( x2_ex ) ];
    r2_rel_v = [ r2_rel_v, norm( b2 - A2 * x2 ) / norm( b2 ) ];
    c2_v = [ c2_v, cond( A2 ) ];
end
figure( 11 ); % for A_1
semilogy( n_vect, e1_rel_v, '-ob', n_vect, r1_rel_v, '-sr', ...
          n_vect, c1_v, '-xk'  );
legend('e_{rel}', 'r_{rel}', 'cond(A_1)');
figure( 12 );  % for A_2
semilogy( n_vect, e2_rel_v, '-ob', n_vect, r2_rel_v, '-sr', ...
          n_vect, c2_v, '-xk'  );
legend('e_{rel}', 'r_{rel}', 'cond(A_2)');