function [x,iter] = mygs(A, b, x0, epsilon, iter_max)%  Gauss-Seidel 迭代解線性方程組%  輸入：A 為系數矩陣， b為右端項，x0為初值，%        epsilon 為相對誤差精度，iter_max 為最大迭代步數%  輸出：x 為解，iter 為迭代步數if prod(diag(A))==0    disp('Error: diag(A) is not invertible!')    returnendn = length(b);x_old = reshape(x0,n,1);x_new = zeros(n,1);for k = 1 : iter_max    for i = 1 : n        tmp = b(i) - A(i,[1:i-1])*x_new([1:i-1]) - A(i,[i+1:n])*x_old([i+1:n]);        x_new(i) = tmp/A(i,i);     end    error1 = norm(x_new - x_old);    if error1 < epsilon, break, end    x_old = x_new;enditer = k;    x = x_new;% END