function y=GS(a,b,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
B=(D-L)\U;
g=(D-L)\b;
y=B*x0+g;
n=1;
while norm(y-x0)>=1.0e-6
    x0=y;
    y=B*x0+g;
    n=n+1;
end
y
n

a=[2 -1 0;-1 3 -1;0 -1 2];
b=[1;8;-5];
GS(a,b,[0;0;0])

y =
    2.0000
    3.0000
   -1.0000
n =
    14