function gradient2
f=@(x) 100*(x(2)-x(1)^2)^2+(x(1)-1)^2;
fp=@(x) [-400*(x(2)-x(1)^2)*x(1)+2*(x(1)-1); 200*(x(2)-x(1)^2)];
[x,xk]=SteepestDescent(f,fp,[0;0]);
F=@(x,y) 100*(y-x.^2).^2+(x-1).^2; % for plotting purposes
Fp1=@(x,y) -400*(y-x.^2).*x+2*(x-1);
Fp2=@(x,y) 200*(y-x.^2);
[X,Y]=meshgrid(-1.5:0.1:1.5,-0.5:0.1:1.5);
contour(X,Y,F(X,Y),50)
hold on
quiver(X,Y,Fp1(X,Y),Fp2(X,Y),5)
plot(xk(1,:),xk(2,:),'-o')
hold off
end




function [x,xk]=SteepestDescent(f,fp,x0,tol,maxiter,tau,be,alinit)
% STEEPESTDESCENT steepest descent min. search with Armijo line search
% [x,xk]=SteepestDescent(f,fp,x0,tol,maxiter,tau,be,alinit) finds an
% approximate minimum of the function f with gradient fp, starting
% at the initial guess x0. The remaining parameters are optional and
% default values are used if they are omited. xk contains all the
% iterates of the method.
if nargin<8, alinit=1; end
if nargin<7, be=0.5; end
if nargin<6, tau=0.99; end
if nargin<5, maxiter=40000; end
if nargin<4, tol=1e-8; end
x=x0;
xk=x0;
p=-feval(fp,x);
k=0;
while norm(p)>tol && k<maxiter
al=alinit;
while feval(f,x+al*p)>feval(f,x)-al*be*p'*p
    al=tau*al;
end
x=x+al*p;
% x=x+0.0009*p;
p=-feval(fp,x);
k=k+1;
xk(:,k+1)=x;
end
niter=k
solution=x
end