function Newton2
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)];
fpp=@(x) [-400*x(2)+1200*x(1)^2+2, -400*x(1); -400*(x(1)), 200];
[x,xk]=Newton(f,fp,fpp,[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,-1:0.1:1.5);
contour(X,Y,F(X,Y),50,LineWidth=1)
hold on
% quiver(X,Y,Fp1(X,Y),Fp2(X,Y),5)
plot(xk(1,:),xk(2,:),'-o')
hold off
end







function [x,xk]=Newton(f,fp,fpp,x0,tol,maxiter,tau,be,alinit)
% NEWTON Minimization with Newton descent and Armijo line search
% [x,xk]=Newton(f,fp,fpp,x0,tol,maxiter,tau,be,alinit) finds an
% approximate minimum of the function f with gradient fp and Hessian
% fpp, 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<9, alinit=1; end
if nargin<8, be=0.1; end
if nargin<7, tau=0.6; end
if nargin<6, maxiter=50000; end
if nargin<5, tol=1e-18; end
x=x0;
xk=x0;
p=-feval(fpp,x)\feval(fp,x);
k=0;
while norm(feval(fp,x))>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+p;
p=-feval(fpp,x)\feval(fp,x);
k=k+1;
xk(:,k+1)=x;
end
k
x
end