Implementation of the Jacobi iteration to find the solution to $Ax = b$

I implemented the Jacobi iteration using Matlab based on this paper, and the code is as follows:

function x = jacobi(A, b)
% Executes iterations of Jacobi's method to solve Ax = b.

assert(size(A, 1) == size(A, 2))
assert(size(A, 1) == size(b, 1))

% matrix with all zeros except the diagonal elements,
% which are those of A.
D = diag(diag(A));

% matrix A but with its diagonal elements equal zero
R = D - A;

tol = 1.e-8;
max_num_of_iter = 1000;

% initialize relative error to large value
relative_error = inf;
num_of_iter = 1;

% initial guess
x_prev = zeros(size(b, 1), 1);

% iterate until convergence or maximum number of iterations is reached
while relative_error > tol && num_of_iter < max_num_of_iter
x = D \ (b + R * x_prev);

relative_error = norm(x - x_prev, inf) / norm(x, inf);

x_prev = x;

num_of_iter = num_of_iter + 1;
end

Now what I don't understand is why the iteration x is found using D \ (b + R * x_prev), more specifically why do we have a + sign instead of a - sign in D \ (b + R * x_prev)?

According to [Wikipedia](https://en.wikipedia.org/wiki/Jacobi_method, the Jacobi method is as follows:

$$\mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)})$$

where $D$ is an $n \times n$ matrix whose entries are all zeros except the diagonal entries which are the diagonal entries of $A$. $R = D - A$, i.e. a matrix whose diagonal elements are zeros and all other are the same as the corresponding elements in $A$.

Clearly, if we want to solve

$$Ax = b$$

and $$A = L + D + U =\\ D + L + U = \\ D + R \Rightarrow \\ A - D = R = \\ L + U$$
then

$$(D + L + U) x = b \\ Dx + Lx + Ux = b \Rightarrow\\ Dx = b - Lx - Ux \Rightarrow \\ Dx = b - (L + U) x$$

since we have defined $L + U = R$, then

$$Dx = b - (L + U) x \Rightarrow \\ Dx = b - Rx \Rightarrow \\ x = D^{-1} \times (b - Rx)$$

So the original formula of the iterative method $\mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)})$ provided by Wikipedia seems to be consistent with my calculations. The problem is that if I use a - instead of + I don't find the right solutions. How do I know that? I'm using x = A \ b (i.e. guassian elimination) to check for the solution.

In your case (line 12 in your code) $$A = D-R.$$ Then \begin{aligned} A\mathbf{x} {}&{}= D\mathbf{x}-R\mathbf{x}\\ {}&{}= \mathbf{b}\,, \end{aligned} which rearranges to $$\mathbf{x} = D^{-1}\left( \mathbf{b} + R\mathbf{x} \right)\,,$$ which leads to the iteration which is the the first line in your while loop.
This only differ's from the Wikipedia version by the definition, i.e. Wikipedia's $$R' = A - D$$ to your $$R = D - A\,.$$