%% Jacobi Iteration

display('Jacobi Iteration');

A=[4,-1,0,1,0;-1,4,-1,0,1;0,-1,4,-1,0;1,0,-1,4,-1;0,1,0,-1,4];
b=[100 100 100 100 100]';

%A=diag(-2.25*ones(7,1)) + diag(ones(6,1),-1) + diag(ones(6,1),1);
%b=[zeros(6,1); -100];

toler=eps;
maxiter=100;
n=size(b);

xk=zeros(n);
x0=ones(n);
k=0;

%Loop until converged, or failed
while ((max(abs((xk-x0)./x0))>toler) & (k<maxiter))
    x0=xk;          %Retain prior solution
    k=k+1;          %Increment counter for iterations
    R=b-A*x0;       %Note - Just a matrix multiply
    xk=x0+R./diag(A);
    display(['k=',int2str(k),'|',num2str(xk','  %1.16e')]);
end;

R
xk