Jacobi Linear System Solver

Jacobi linear system solver. For more info, consult the internet: http://en.wikipedia.org/wiki/Jacobi_method

function [k,x] = Jacobi(n,A,b)
% Jacobi linear system solver
% Where A is a square matrix of size n by n and b is a 1xn matrix

k_max = 100;
delta = 1.0e-10;
eps = 1.0e-5;
x=zeros(n,1);
for k=1:k_max
    y = x;
    for i = 1:n               
        sum = b(i);
        diag = A(i,i);
        if abs(diag) < delta
            error('diagonal element too small')
        end
        
        for j = 1:n 
            if j ~= i
                sum = sum - A(i,j)*y(j);
            end
        end
        
        x(i) = sum/diag;
    end
    
    err = x - y;
    err=abs(err);
    if ( max(abs(err)) < eps ) %use infinity norm 
        break; % convergence achieved
    end
end

if k == k_max
    disp('maximum iterations reached')
end

% end Jacobi