% m22_pureNewtonmin.m - minimisation of f(x,y) by pure Newton method

 
  f = @(x,y) x.^2+y.^2+y.*sin(x)-sin(2*y);    % function

  g = @(x,y) [      2*x+y.*cos(x)             % gradient
               2*y+sin(x)-2*cos(2*y) ];

  H = @(x,y)  [2-y.*sin(x)      cos(x)        % Hessian
                 cos(x)      2+4*sin(2*y)];

  x = [-3:.1:3]; y = x;
  [xx,yy] = meshgrid(x,y);               % set up plotting grid
  ff = f(xx,yy);                         % evaluate function on grid
  hold off, contour(xx,yy,ff,15)         % make contour plot
  colorbar, grid on, hold on
  LW = 'linewidth'; MS = 'markersize';
  plot(-.2714,.5635,'.k',MS,30)          % the solution: big dot
  x = input('initial x? ');
  y = input('initial y? ');
  xold = x; yold = y;
  for k = 0:15
    plot([xold x],[yold y],'-r',LW,1.2)  % plot line to current point
    plot(x,y,'.r',MS,15)                 % plot current point
    disp(' '), pause
    fxy = f(x,y)
    gxy = g(x,y)                         % evaluate f, g, H
    Hxy = H(x,y);
    xold = x; yold = y;
    s = -Hxy\gxy;                        % compute next step
    x = x+s(1); y = y+s(2);              % update x and y
  end