   function [x,nrmHist] = steep(A, b, x, nsteps, tol) 
%  function [x,nrmHist] = steep(A, b, x, nsteps, tol) 
%  steepest descent algorithm for solving Ax=b
%  A = matrix  [must be SPD] 
%  x = initial guess on entry -- app sol on return 
%  b = rhs, 
%  nsteps = max number of steps --
 r = b - A * x;
 ro = norm(r)
 nrmHist(1) = ro;
 if (nargin<5), tol = 1.e-08;, end
 tol_rel = tol*ro;
%%-------------------- loop
 for i=1:nsteps
   ar = A * r; 
   alp = r'*r / ( r' * ar);
   x = x + alp * r;
   r = r - alp* ar;
   ro = norm(r);
%%--------------------  save norms for plotting
   nrmHist(i+1) = ro;
   if (ro < tol_rel), break, end 
 end 
          
