function [GSx, Jx, elteres] = Gauss_Jacobi_it(A,b,maxiter)

S = tril(A);
T = S - A;
GSx = zeros(size(b));
for i = 1:maxiter
	GSx= S\(T*GSx+b);
end

D = diag(diag(A));
T = D-A;
Jx = zeros(size(b));
for k = 1:maxiter
	Jx = D\(T*Jx+b);
end
elteres=norm(GSx-Jx);