function x = Jacobi( A, b, error_limit )

s = size(A,1);
D = diag(diag(A));
L = tril(A) - D;
U = triu(A) - D;
D_inv = inv(D);
B = - D_inv * (L + U);
c = D_inv * b;

x1 = zeros(s,1);
x2 = rand(s,1);
while sum((x2 - x1).^2) > error_limit^2
    x1 = x2;
    x2 = B * x1 + c;
end

x = x2;

end

