
tol = 1e-6;

x0 = [10 10]';

[x, k, nf, fval] = mynlsteepest(@rosen, @rosengrad, x0, tol) 
  
function [x, k, nf, obj] = mynlsteepest(func, grad, x, tol)
%MYNLSTEEPEST implements steepest descent method
%             using Armijo stepsize rule. It terminates 
%             when the norm of the gradient is below tol
  
  % Armijo stepsize rule parameters
  sigma = .1;
  beta = .5;
  
  obj = func(x);
  g = grad(x);
  
  k = 0;   % k = # iterations
  nf = 1;  % nf = # function eval

  while norm(g) > tol
    d = -g;         % steepest descent direction
    a = 1;
    newobj = func(x + a*d);
    nf = nf + 1;
    while (newobj - obj)/a > sigma*g'*d
      a = a*beta;
      newobj = func(x + a*d);
      nf = nf + 1;
    end
    x = x + a*d;
    obj = newobj;
    g = grad(x);
    k = k + 1;
  end
end


function y = rosen(x)
%This routine evaluates the Rosenbrock function.
  y = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end

function y = rosengrad(x)
%This routine evaluates the gradient of the Rosenbrock function
  y(1) = 100*(2*(x(1)^2-x(2))*2*x(1)) + 2*(x(1)-1);
  y(2) = 100*(-2*(x(1)^2-x(2)));
  y = y';
end



