function X=jakobi(A,b,x0,maxiter)
% Jacobi method - solves Ax=b with the given initial approximation x0
% Computes fixed number of iterations - maxiter
% Returns matrix X with columns containing all intermediate solutions

x_k=x0; %initialize
X=x_k; % store the first solution - initial guess 
D=diag(diag(A)); %create diagonal matrix D from the diagonal of A

for i=1:maxiter

    y=D\(b-A*x_k); % compute update
    x_k=x_k+y; % update x_k
    X=[X,x_k]; % store t

end

% Note that this implementation might not be optimal in terms of storage
%  for example storing the diagonal of A to a diagonal matrix D is
%  redundant, but its done this way so that its easier to read the code :)

end