using LinearAlgebra

function jacobi(A,b,x, tol=1E-5, maxit=10000)
    
    for k in 1:maxit
        xold = copy(x)
        R = b - A*xold
        x = xold + R./diag(A)
        relerr = norm(R)/norm(x)
        if relerr <= tol
            return x, k
        end
    end
    println("Warning: no convergence in $maxit iterations")    
    return x, maxit
end;

using Statistics
using PyFormattedStrings

n = 60

A  = diagm(-3=>fill(1,n-3), -1=>fill(-1,n-1), 0=>fill(4,n), 1=>fill(-1,n-1), 3=>fill(1,n-3))
x0 = zeros(n)
b  = ones(n)

x, k_iter = jacobi(A,b,x0)

print(f"The problem required {k_iter} iterations")

x_ju = A\b
reAvg = mean(abs.((x_ju-x)./x_ju))
println(f"The average relative error between the Jacobi function and direct solver is {reAvg:.3e}")
println("\nx =$x")