function [x, flag_converge, iter] = Jacobi(A,b,tol,n_iter,x0)
% x = Jacobi(A,b,tol,n_iter,x0) tries to solve the equation Ax=b using
% Jacobi iteration.
%
% A ... n times n matrix
% b ... n times 1 vector
% tol ... tolerance level, optional, default value = 10^(-6)
% n_iter ... maximal number of iterations, optional, default value = 100
% x0 ... starting vector, optional, default value = 0
%
% additional output arguments:
% flag_converge ... a flag that is 1 if the method achieved the required
%       tolerance and is 0 else
% iter ... actual number of iterations


% set the default starting value to 0 if undefined
if(nargin < 5)
    x0 = zeros(size(b));
end

% set the default number of iterations to 100 if undefined
if(nargin < 4 || isempty(n_iter))
    n_iter = 100;
end

% set the default accuracy to 10^(-6) if underfind
if (nargin < 3 || isempty(tol))
    tol = 10^(-6);
end

% Initialize the result
x = x0;

% define the matrix/vector driving the iteration
B = eye(length(b))-diag(1./diag(A))*A;
c = b./diag(A);

% intialize the convergence flag
flag_converge = 0;

% main iteration
for k=1:n_iter
    % update of x
    x_new = B*x+c;
    % test for convergence
    if(norm(x_new-x) < tol)
        x = x_new;
        flag_converge = 1;
        break;
    end
    % update the iterate
    x = x_new;
end

% store the number of iterations
iter = k;