   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;
 nrmHist(1) = norm(r); 
 if (nargin<5), tol = 1.e-08;, end
 tol1 = tol*nrmHist(1);
%%-------------------- 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 < tol1), break, end 
 end 
          
