function [out, Tgs,cgs] = gauss_seidel(L,U,D,b,x0,n)
x = x0;
Tgs = (D-L)\U;
cgs = (D-L)\b;
for i=2:n
    x(:,i) = Tgs*x(:,i-1) + cgs;
end
out = x;
return;