function[x_k,error_vec]=jacobi(x_k,A,b,max_iter)

n = sqrt(length(x_k));
error_vec = [];
x_true = A\b;
nx_true = norm(x_true);
max_x_true = max(x_true(:));

% Perform Jacobi Iterations to approximately solve Ax=b.

D = diag(diag(A));
L = -tril(A,-1);
U = -triu(A,+1);

disp(' Jacobi Iterations ')
for i = 1:max_iter    
    h = waitbar(i/max_iter);
    % Save relative errors. 
    error_vec = [error_vec, norm(x_k(:)-x_true(:))/nx_true];
    % Plot solution. 
    if mod(i,10) == 0
      subplot(221)
        imagesc(reshape(x_k,n,n))
        title(['Plot of ',num2str(i),' Jacobi iterate'])
        subplot(224)
        imagesc(reshape(x_true,n,n))
        title('Exact Solution')
        drawnow, pause(0.1)
    end
    % Jacobi iteration
    x_k = D\(b+(L+U)*x_k); 
end
close(h)