% compute and plot models of a function R^2 -> R 
%
% To get help for a command "command_name" type
%   help command_name
%
% To get an overview over all help topics type
%   help 
% without any arguments.

%
% specify grid (domain) of evaluation
%

x=0.001:.05:2;      % grid points along x-axes         
y=0.001:.05:2;      % grid points along y-axes

nx=length(x);       % size of grid along x-axes 
ny=length(y);       % siye of grid along y-axes

%
% evaluate function in grid points
%

fz=zeros(nx,ny);    % intiliaze matrix size to grid size and set to zero  
for i=1:nx
  for j=1:ny
    % compute f(x(i),y(j)) and store in fz(i,j)
    fz(i,j)= -log(x(i))-5*(x(i)-.5).^3+7*x(i)-log(y(j))-5*(y(j)-.5).^3+7*y(j);
  end
end

%
% fix a point [xp;yp] 
% its first and second order information gives rise to the model around [xp;yp]
%

xp=x(10);   % choose x-coordinate as a grid point for convenience
yp=y(10);   % choose y-coordinate as a grid point for convenience

fzp=-log(xp)-5*(xp-.5).^3+7*xp-log(yp)-5*(yp-.5).^3+7*yp;  % f(xp,yp)
gradf=[-1./xp-15*(xp-.5).^2+7;-1./yp-15*(yp-.5).^2+7];     % gradient(xp,yp)
Hf=[1./(xp.^2)-30*(xp-.5), 0;0 1./(yp.^2)-30*(yp-.5)];     % Hessian(xp,yp)

m1z=zeros(nx,ny);    % initialize first order model matrix for grid size
m2z=zeros(nx,ny);    % initialize second order model matrix for grid size
for i=1:nx
  for j=1:ny
    dvec=[x(i)-xp;y(j)-yp];                     % point of evaluation
    m1z(i,j)= fzp+gradf'*dvec;                  % value of first order model
    m2z(i,j)= fzp+gradf'*dvec+dvec'*Hf*dvec/2;  % value of second order model
  end
end

%
% plot function, first and second order model
%

figure(1);               % output next plot on window of figure 1 
hold off;                % overwrite old plot in figure 1
color0=zeros(nx,ny);     % set color in all grid points to color of value 0
mesh(x,y,fz,color0);     % plot the function
hold on;                 % keep old plot while adding further plots
color1=ones(nx,ny);      % set color of next plot to color value 1 
mesh(x,y,m1z,color1);    % plot the first order model
color2=2*ones(nx,ny);    % set color of next plot to color value 2
mesh(x,y,m2z,color2);    % plot the second order model
view(70,30);             % set the view angle for the plot
axis([0 2 0 2 -10 10]);  % set value range to plot
%print some title for the plot giving the coordinates [xp;yp] 
title([' function, first and second order model at x=',num2str(xp),', y=',num2str(yp)]);

