A=[4 -1 0; -1 4 -1; 0 -1 4];       
M=diag(diag(A));                     % decomposition A=M-N
N=M-A;   
G=M\N;                               % iteration matrix
rho=max(abs(eig(G)))                 % spectral radius of G
R_inf=-log(rho)                      % asymptoric convergence rate
delta=1000;                          % reduction factor
epsi=1/delta;                        % tolerance
K=-log(epsi)/R_inf                   % estimate of the iterations "k"
e=ones(3,1);                         % initial error
nor(1)=norm(e,2);
for j=1:round(K)                     % approximate K by the closest integer
  e=G*e;                             % iteration step
  nor(j+1)=norm(e,2)/nor(1);
  disp(['j=' num2str(j) ' | norm=' num2str(nor(j+1))]);
end
v=0:round(K);
semilogy(v,nor,'ob-',v,10*rho.^v,'b');
legend('error','asymptotic convergence');
xlabel('iteration');
ylabel('error'); 
set(gca,'fontsize',16);                   % bigger font size
