function z = gs1(A,b,x0,N)
% GS1 Do N steps of Gauss-Seidel iteration. Implemented using
% Matlab tools including backslash.
% Usage:
% >> z = gs1(A,b,x0,N)
DL = tril(A); % A = (D+L) + U; DL includes diagonal
U = triu(A,1); % super-diagonal
z = x0;
for k = 1:N
z = DL \ (b - U * z);
end