%  Algorithm 1.  Simple Iteration.

%  User enters A, M, b, x0, and itmax.
%  Prints out residual norm at each step and also
%  returns residual norms in array resid.

n = length(b);
resid = zeros(itmax+1,1);

xk = x0;
rk = b - A*xk;
rknrm = norm(rk),
resid(1) = rknrm;
zk = M\rk;

for k=1:itmax,
  xk = xk + zk;
  rk = b - A*xk;
  rknrm = norm(rk),
  resid(k+1) = rknrm;
  zk = M\rk;
end;