n1 = 3;
A1 = [ 3 -2 1; 2 1.65 -1; 0 1 4 ];
P1_J = diag( diag( A1 ) );  % preconditioning matrix Jacobi
B1_J = eye( n1, n1 ) - inv( P1_J ) * A1; % iteration matrix Jacobi 
rho1_J = max(abs(eig( B1_J ))),   % spectral radius B1_J
% rho1_J =
%     0.9851
P1_GS = diag( diag( A1 ) ) - ( - 1 * tril( A1, -1 ) ); % prec.matr.Gauss-Seidel
B1_GS = eye( n1, n1 ) - inv( P1_GS ) * A1; % iteration matrix Gauss-Seidel
rho1_GS = max(abs(eig( B1_GS ))),   % spectral radius B1_GS
% rho1_GS =
%     1.0606

n2 = 3;
A2 = [ 5 -3 -2; -3 3 0; -2 0 4 ];
eig_A2 = eig(A2)'
% eig_A2 =
%     0.4103    3.7126    7.8771

P2_J = diag( diag( A2 ) ); % Jacobi
B2_J = eye( n2, n2 ) - inv( P2_J ) * A2;
rho2_J = max(abs(eig( B2_J )))
% rho2_J =
%     0.8944
P2_GS = diag( diag( A2 ) ) - ( - 1 * tril( A2, -1 ) ); % Gauss-Seidel
B2_GS = eye( n2, n2 ) - inv( P2_GS ) * A2;
rho2_GS = max(abs(eig( B2_GS )))
% rho2_GS =
%     0.8000

n3 = 100;
A3 = diag( 4 * ones( n3, 1 ), 0 ) + diag( - 1 * ones( n3 - 1, 1 ), 1 ) + ...
     diag( - 1 * ones( n3 - 1, 1 ), -1 );
P3_J = diag( diag( A3 ) );  % Jacobi
B3_J = eye( n3, n3 ) - inv( P3_J ) * A3;
rho3_J = max(abs(eig( B3_J )))
% rho3_J =
%     0.4998
P3_GS = diag( diag( A3 ) ) - ( - 1 * tril( A3, -1 ) ); % Gauss-Seidel
B3_GS = eye( n3, n3 ) - inv( P3_GS ) * A3;
rho3_GS = max(abs(eig( B3_GS )))
% rho3_GS =
%     0.2498

x1 = ones( n1, 1 );     b1 = A1 * x1;    x0 = zeros( n1, 1 );
[ x1_J, k1_J, res1_J ] = jacobi( A1, b1, x0, 1e-6, 1000 );
err1_J = norm( x1 - x1_J ),    k1_J,    res1_J
% err1_J =
%    2.7745e-07
% k1_J =
%    969
% res1_J =
%    9.6699e-07

x2 = ones( n2, 1 );   b2 = A2 * x2;   x0 = zeros( n2, 1 );
[ x2_J, k2_J, res2_J ] = jacobi( A2, b2, x0, 1e-6, 1000 );  % Jacobi
err2_J = norm( x2 - x2_J ),    k2_J,   res2_J
% err2_J =
%    1.6855e-06
% k2_J =
%    124
% res2_J =
%    8.8408e-07
[ x2_GS, k2_GS, res2_GS ] = gauss_seidel( A2, b2, x0, 1e-6, 1000 ); % Gauss-S.
err2_GS = norm( x2 - x2_GS ),   k2_GS,   res2_GS
% err2_GS =
%    1.4712e-06
% k2_GS =
%    63
% res2_GS =
%    9.8080e-07

x3 = ones( n3, 1 );    b3 = A3 * x3;    x0 = zeros( n3, 1 );
[ x3_J, k3_J, res3_J ] = jacobi( A3, b3, x0, 1e-6, 1000 );  % Jacobi
err3_J = norm( x3 - x3_J ), k3_J, res3_J
% err3_J =
%    2.8220e-07
% k3_J =
%    25
% res3_J =
%    5.6510e-07
[ x3_GS, k3_GS, res3_GS ] = gauss_seidel( A3, b3, x0, 1e-6, 1000 ); % Gauss-S.
err3_GS = norm( x3 - x3_GS ), k3_GS, res3_GS
% err3_GS =
%    2.1655e-07
% k3_GS =
%    16
% res3_GS =
%    4.3476e-07

gamma_v = linspace( -10, 25, 2001 );
rho4_J_v = [ ];   rho4_GS_v = [ ]; 
for gamma = gamma_v 
    n4 = 4;
    A4 = [ 8 gamma -2 -1; -2 2 -gamma -3; -1 -2 18 -18; -1 -3 -7 25 ];
    P4_J = diag( diag( A4 ) );  % Jacobi
    B4_J = eye( n4, n4 ) - inv( P4_J ) * A4;
    rho4_J_v = [ rho4_J_v, max( abs( eig( B4_J ) ) ) ];
    P4_GS = diag( diag( A4 ) ) - ( - 1 * tril( A4, -1 ) );  % Gauss-S.
    B4_GS = eye( n4, n4 ) - inv( P4_GS ) * A4;
    rho4_GS_v = [ rho4_GS_v, max( abs( eig( B4_GS ) ) ) ];
end
plot( gamma_v, rho4_J_v, '-b', gamma_v, rho4_GS_v, '-r', ...
    gamma_v, 1 + 0 * gamma_v, '--k' );
legend('\rho Jacobi','\rho Gauss-Seidel','1' );