% Jacobi scheme implementation using matrix splitting A=P-Q

A=[7 2 -3 1; -4 8 1 -2; 3 -2 9 2; 1 -1 1 4];
b=[-1;16;-3;5];

x=[1;1;1;1]; % initial guess

itermax=50;
tol=1e-8;

n=length(b);

% set up P and Q
P=diag(diag(A));
Q=diag(diag(A))-A;

for k=1:itermax
    % assign most recent x to xprev
    xprev=x;
    % calculate new x by solving P*x=Q*xprev+b by backward or forward substitution
    x=backsub(P,Q*xprev+b);   
    % print iteration number and the x just computed
    fprintf('%2d',k)
    for i=1:n-1
        fprintf('%12.8f  ',x(i))
    end
    fprintf('%12.8f\n',x(n))
    % check convergence
    if norm( x - xprev ) <= tol
        break
    end
end