function gauss_seidel(a, b, iterations, x = zeros(length(b)))
# The matrix should not have near-0 values in the diagonal
# (Ideally it should be diagonally dominant for certain convergence)
n = length(x)
for i = 1:iterations
for j = 1:n
x[j] = b[j]
for k = 1:n
if j != k
x[j] -= a[j,k] * x[k]
end
end
x[j] /= a[j,j]
end
end
x
end