function [x, B, K, residuo] = metodi_iter(A,b)

% METODI ITERATIVI
% OUTPUT 
% x: soluzione appross. del sistema Ax=b
% B: matrice di iterazione
% K: numero di iterazioni
% residuo: ||Ax-b||

disp('Quale metodo iterativo vuoi utilizzare?')
flag=input('Jacobi --> digita J, Gauss Seidel --> digita GS ','s') ;
n=size(A,2) ;
x=zeros(n,1) ;
invD=eye(n) ;
K=0 ;
errore=10 ;  
normb=norm(b) ;    
tol=input('Inserire il valore massimo dell''errore ' ) ;
maxiter=input('Inserire il numero massimo di iterazioni: ') ;
err=zeros(maxiter,1) ;  % err=[] ; vettore degli errori

if flag=='J'
    for i=1:n
        if A(i,i)~=0
          invD(i,i)=1/A(i,i) ;
        else
          disp('Jacobi non applicabile')
          return
        end
    end
    B=eye(n)-invD*A ;
elseif strcmp(flag,'GS')
    M=tril(A) ;
    B=eye(n)-(M\eye(size(M)))*A ;
end
    eig(B)
    rho=max(abs(eig(B))) ;
    disp(['RAGGIO SPETTRALE DI B: ', num2str(rho)])

tic
% CALCOLO DELL'ITERATA x_(k+1) e test d'arresto
while errore>tol && K<=maxiter
    if flag=='J'
        xold=x ;
        for i=1:n
            x(i)=(b(i)-A(i,1:i-1)*xold(1:i-1)-A(i,i+1:n)*xold(i+1:n))/A(i,i) ;
        end
    elseif strcmp(flag,'GS')
         for i=1:n
             x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x(i+1:n))/A(i,i) ;
         end
    else
        B=zeros(n) ;
        return
    end
    residuo=A*x-b ;
    errore=norm(residuo)/normb ; %stimatore dell'errore relativo al passo K
    K=K+1 ;
    err(K)=errore ;
end
t=toc ;
disp(['Tempo di calcolo: ', num2str(t)])
disp(['Numero di iterazioni effettuate: ' num2str(K)])
residuo=norm(A*x-b) ;
figure
semilogy(err(1:K),'--')
grid
title(['Andamento dell''errore con il metodo di ', flag , ' - Errore finale=' num2str(errore)])
end