Jacobi iteration code, not producing correct solution

I have written a code for Jacobi iteration but the solution coming through my code and matlab's A\b is not giving me the same solution. I do not find any mistake in my code. Any help appreciated.

function f=Jacobi(A,b)
%Breaking the matrix in diagonal, Triangular and upper triangular
D=diag(diag(A));
L=tril(A,-1);
U=triu(A,1);

%Forming the Jacobi matrix
M=D;N=-(L+U);
J=inv(M)*N;
c=inv(M)*b;

%specifying the error Bound
Err=1e-6;

%doing the Jacobi iteration
n=length(A);
x=[0 -1 1]';
x=J*x+c;
x0=x;
while norm(x0-x)>Err
  x0=x;
  x=J*x+c;
end
f=x;
end

Answer 1

[David Goodmanson]

Hi Arindam,

Just before the while loop you need to reverse the two statements to make it

 x0=x;
 x=J*x+c;

With the order you have them, the while statement is satisfied automatically. Also your initial x vector is hardwired right now to length 3, better would be something like x = rand(n,1).