function [U,c] = equi_system(A,b)

v = diag(A); % Vettore con gli elementi diagonali di A

if (det(A)==0) | (v(find(v))~=v) | (size(A,1)~=size(A,2)) | size(A,1)~=size(b,1) % Utilizzo di un OR vettoriale
    error('La matrice A o il vettore b non soddisfano le condizioni per applicare il metodo')
end

n = size(A,1); % Dimensione di A (quadrata)

% Inizializzo le variabili
Akp = zeros(size(A));
Ak = A;
c = b;

for i= 2:n % Ciclo sulle iterate
    Akp(i-1,i-1:n) = Ak(i-1,i-1:n); % Riga i-1 non viene cambiata
    Akp(i:n,i-1) = 0; % Colonna sotto l'elemento i-1 viene posta a zero
    m = Ak(i:n,i-1)/Ak(i-1,i-1);  % Aggiorno gli altri elementi della matrice
    Akp(i:n,i:n) = Ak(i:n,i:n) - m.*Ak(i-1,i:n);
    
    % Aggiorno il termine noto
    c(i:n) = c(i:n) - m*c(i-1);
    % Pongo la matrice nuova come quella da cui partire nella nuova iterata
    Ak = Akp;
end

% Definisco U, output della funzione
U = Akp;


end