function s=Gauss(a,b,x0,eps,N,)
% Gauss-Seidel迭代法解線性方程組
% a為系數矩陣，b為方程組ax=b的右端值，x0為初值
if nargin==4
    eps=1.0e-6;
elseif nargin<4
%    error
    return
end
D=diag(diag(a));%求對角矩陣
L=tril(a,-1);%求嚴格下三角矩陣
U=triu(a,1);%求嚴格上三角矩陣
C=inv(D+L);
B=-C*U;
f=C*b;
s=B*x0+f;
i=0;
while (norm(s-x0)>=eps&&(i<=N))
    x0=s;
    s=B*x0+f;
    i=i+1;
end
return