APL(NARS), 101 chars

r←a f w;b;x;e;D;R;A;t
(b x e)←w⋄R←a-D←a×=/¨⍳⍴a⋄A←⌹D⋄t←b⌹a⋄r←0
→3×⍳e>⌈/∣x-t⋄x←A+.×b-R+.×x⋄r+←1⋄→2
r,←⊂x

//22+39+35+5=101 it would return the number of iteration + the solution found.
For me because there is a convergience in the result, the way to go would be find the difference from this computation x and the precedent one, and calculate the max of abs of the differences , if that number is <e than stop and return the last x. This should be ok only if the ⌈/|x_i-x_(i+1) is decreasing toward 0 in the increase of i iteration.

Test

  (2 2 ⍴ 9 ¯2 1 3)f (3 4)  (1 1) 0.04
2  0.5555555556 1.148148148 
  (2 2 ⍴ 2 3  1 4) f (2 ¯1) (2.7 ¯0.7)  1
0  2.7 ¯0.7 
  (2 2 ⍴ 9 ¯2 1 3)f (3 4)  (1 1) 0.000004
10  0.5862059737 1.137931342 
  

APL (Dyalog), 78 68 65 49 bytes

Exactly the type of problem APL was created for.

-3 thanks to Erik the Outgolfer. -11 thanks to ngn.

Anonymous infix function. Takes A b e as left argument and x as right argument. Prints result to STDOUT as vertical unary using 1 as tally marks, followed by 0 as punctuation. This means that even a 0-result can be seen, being no 1s before the 0.

{⎕←∨/e<|⍵-b⌹⊃A b e←⍺:⍺∇D+.×b-⍵+.×⍨A-⌹D←⌹A×=/¨⍳⍴A}

Try it online!

Explanation in reading order

Notice how the code reads very similarly to the problem specification:

{…} on the given A, b, and e, and and the given x,
 ⎕← print
 ∨/ whether there is any truth in the statement that
 e< e is smaller than
 |⍵- the absolute value of x minus
 b⌹ b matrix-divided by
 ⊃A b e the first of A, b, and e (i.e. A)
 ←⍺ which are the left argument
 : and if so,
  ⍺∇ recurse on
  D+.× D matrix-times
  b- b minus
  ⍵+.×⍨ x, matrix multiplied by
  A- A minus
  ⌹D the inverse of D (the remaining entries)
  ← where D is
  A× A where
  =/¨ there are equal
  ⍳ coordinates for
  ⍴A the shape of A (i.e. the diagonal)

Step-by-step explanation

The actual order of execution right-to-left:

{…} anonymous function where ⍺ is A b e and ⍵ is x:
 A b c←⍺ split left argument into A, b, and e
 ⊃ pick the first (A)
 b⌹ matrix division with b (gives true value of x)
 ⍵- differences between current values of x and those
 | absolute values
 e< acceptable error less than those?
 ∨/ true for any? (lit. OR reduction)
 ⎕← print that Boolean to STDOUT
 : and if so:
  ⍴A shape of A
  ⍳ matrix of that shape where each cell is its own coordinates
  =/¨ for each cell, is the vertical and horizontal coordinates equal? (diagonal)
  A× multiply the cells of A with the that (extracts diagonal)
  ⌹ matrix inverse
  D← store in D (for Diagonal)
  ⌹ inverse (back to normal)
  A- subtract from A
  ⍵+.×⍨ matrix multiply (same thing as dot product, hence the .) that with x
  b- subtract that from b
  D+.× matrix product of D and that
  ⍺∇ apply this function with given A b e and that as new value of x

Clojure, 212 198 196 bytes

#(let[E(range)I(iterate(fn[X](map(fn[r b d](/(- b(apply +(map * r X)))d))(map assoc % E(repeat 0))%2(map nth % E)))%3)](count(for[x I :while(not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))]x)))

Implemented without a matrix library, it iterates the process 1e9 times to get the correct answer. This wouldn't work on too ill-conditioned inputs but should work fine in practice.

Less golfed, I was quite happy with the expressions of R and D :) The first input % (A) has to be a vector, not a list so that assoc can be used.

(def f #(let[R(map assoc %(range)(repeat 0))
             D(map nth %(range))
             I(iterate(fn[X](map(fn[r x b d](/(- b(apply +(map * r x)))d))R(repeat X)%2 D))%3)]
          (->> I
               (take-while (fn[x](not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))))
               count)))

R, 138 bytes

function(A,x,b,e){R=A-(D=diag(diag(A)))
g=solve(A,b)
if(norm(t(x-g),"M")<e)T=0
while(norm((y=solve(D)%*%(b-R%*%x))-g,"M")>e){T=T+1;x=y}
T}

Try it online!

thanks to NikoNyrh for fixing a bug

It's also worth noting that there is an R package, Rlinsolve that contains a lsolve.jacobi function, returning a list with x (the solution) and iter (the number of iterations required), but I'm not sure that it does the correct computations.

Python 3, 132 bytes

f=lambda A,b,x,e:e<l.norm(x-dot(l.inv(A),b))and 1+f(A,b,dot(l.inv(d(d(A))),b-dot(A-d(d(A)),x)),e)
from numpy import*
l=linalg
d=diag

Try it online!

Uses a recursive solution.