n = 10; 
A = 3*diag(ones(n,1)) - ...
    (diag(ones(n-1,1),-1)+diag(ones(n-1,1),1));
b = sum(A,2);
x = ones(n,1);  

% A = L + D + U

L = tril(A,-1);
D = diag(diag(A));
U = triu(A,1);

% niters 

niters = 25;

% Jacobi

MJ = -inv(D)*(L+U);
bJ =  inv(D)*b;
xJ  = zeros(n,niters+1);

% GS

MGS   = -inv(L+D)*U;
bGS   =  inv(L+D)*b;
xGS   = zeros(n,niters+1);


% compute iterations

for k = 1:niters
    xJ(:,k+1)   = MJ*xJ(:,k) + bJ;
    xGS(:,k+1)  = MGS*xGS(:,k) + bGS;
end

%
% display error & error ratios & spec radius
%

format short e

% 2-norm

errJ    = vecnorm(repmat(x,1,niters+1)-xJ)';
errGS   = vecnorm(repmat(x,1,niters+1)-xGS)';

disp(' ')
disp('    err(MJ)      err(MGS)  ')
disp('   ----------   ----------')

disp([errJ errGS])
disp(' ')

disp('   eratio(MJ)   eratio(MGS) ')
disp('   ----------   ---------- ')
disp([errJ(2:niters+1)./errJ(1:niters) ...
      errGS(2:niters+1)./errGS(1:niters)]);

disp(' ')
disp('    rho(MJ)      rho(MGS)  ')
disp('   ----------   ----------')
rhoJ   = max(abs(eig(MJ)));
rhoGS  = max(abs(eig(MGS)));
disp([rhoJ rhoGS])


